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Abstract 

In this paper a family of fixed point algorithms, generalizing the Petviashvili method, 
is considered. A previous work studied the convergence of the methods. Presented 
here is a second part of the analysis, concerning the introduction of some accelera¬ 
tion techniques into the iterative procedures. The purpose of the research is two-fold: 
one is improving the performance of the methods in case of convergence and the sec¬ 
ond one is widening their application when generating traveling waves in nonlinear 
dispersive wave equations, transforming some divergent into convergent cases. Two 
families of acceleration techniques are considered: the vector extrapolation meth¬ 
ods and the Anderson acceleration methods. A comparative study through several 
numerical experiments is carried out. 
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1 Introduction 


In a previous paper [3], a family of fixed-point algorithms for the numerical 
approximation of nonlinear systems of the form 


Lu = N{u), u G M™', m > 1, (1) 

was introduced. In ([1]), L is a nonsingular mxm real matrix and N : —)■ M”* 

is an homogeneous function with degree p, \p\ > 1 (this means that N{Xu) = 
XPN{u)). Among other applications, systems of this form are very typical in 
the numerical generation of traveling waves and ground states in water wave 
problems and nonlinear optics. For the numerical approximation to solutions 
of (11]), the use of the classical fixed-point algorithm is not suitable. This is due 
to the fact that if u* is a solution and S = L~^N'{u*) stands for the iteration 
matrix at u*, then, since N is homogeneous of degree p then N'{u*)u* = 
pN{u*) and therefore 


Su* = L-^N\u*)u* = pL-^N{u*) = pu*, 


that is, p is an eigenvalue of S with magnitude above one. This makes the 
iteration not convergent in general. 


As an alternative and based on the Petviashvili method, [56115 5 1146 IIT7] . the 
following fixed-point algorithms were considered in [3]: 


LUn+l = s{Un)N{Un), n = 0, 1, . . . , (2) 

from mq 7 ^ 0 and where s : —)■ M is a function satisfying the following 

properties: 

(PI) A set of fixed points of the iteration operator 

F{x) = s{x)L-^N{x), (3) 

coincides with a set of fixed points of ([T|). This means that: (a) if u* is a 
solution of ([T]) then s{u*) = 1; (b) inversely, if the sequence {Mn}n, generated 
by ([2|), converges to some ?/, then s{y) = 1 (and, consequently, y is a. solution 
of (HD). 

(P2) s is homogeneous with degree q such that \p + q\ < 1. 

The function s is called the stabilizing factor of the method, inheriting the 
nomenclature of the Petviashvili method. Actually, formula (ED generalizes the 
Petviashvili scheme, which corresponds to the choice 
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{Lx, x) 


7 


{N{x),x) 


, 9 = 7(1-P), 


( 4 ) 


The first part of the work, carried out in [3] (see also 0 ). analyzed the conver¬ 
gence of (E]). The main conclusion was that, compared to the classical hxed- 
point algorithm, the stabilizing factor acts like a hlter of the spectrum of the 
matrix S in the sense that: 

• The eigenvalue A = p of S' is transformed to the eigenvalue A = p -|- g of the 
iteration matrix F'[u*) of ([3]). 

• The rest of the spectrum of F'{u*) is included into the spectrum of S. 

Thus, the convergence of ([2]) depends of the spectrum of S different from p. 
From these conclusions, several results of convergence can be derived, see [3] 
for details. 

The motivation of this paper is two-fold. First, several numerical experiments 
in the literature show that the Petviashvili type algorithms are in sometimes 
computationally slower than other alternatives. In order to continue to bene¬ 
fit from the easy implementation (one of the advantages of the methods) the 
algorithms should improve their performance with the inclusion of some ac¬ 
celeration technique. A further motivation comes from the known mechanism 
of some extrapolation methods, |6l], to transform divergent into convergent 
cases. The application of this property to these Petviashvili type methods 
may extend their use to compute traveling waves under more demanding con¬ 
ditions, for example in two dimensions or/and in case of highly oscillatory 
waves. 

The literature on acceleration techniques is very rich with many different fam¬ 
ilies and strategies, IMIEZ] . This paper will be focused on two types of pro¬ 
cedures: the vector extrapolation methods, [2^70ll41ll20ll66j and the Anderson 
mixing, PI75II53II34] . We think that the first one is the most widely studied 
group; in particular, known convergence results for some of these methods 
will serve us to justify several examples of transformation from divergence 
to convergence when generating traveling waves iteratively. The second family 
accelerates the convergence by introducing the strategy of minimization of the 
residual in some norms at each step. They have been revealed efficient in, for 
example, electronic structure computations, pf58] (see also [531175] and refer¬ 
ences therein) and, to our knowledge, this is the hrst time they are applied to 
the numerical generation of traveling waves. 

The main purpose of this paper is then exploring by numerical means the 
application of these two acceleration methods to the generation of traveling 
waves through several problems of interest and from the Petviashvili type 
methods ([2]) and their extended versions derived in [1|. With the case studies 
presented here we have tried to cover different situations of hard computation 
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of the waves as an attempt to give some guidelines of application. In this sense, 
the paper provides several conclusions to be emphasized: 

• The use of acceleration techniques is highly recommended here since it im¬ 
proves the performance in general and allows to extend the application of 
the methods to computationally harder situations, with especial emphasis 
on two-dimensional simulations and highly oscillatory wave generation. 

• By comparing the two families of acceleration techniques considered in this 
study, the vector extrapolation methods are in general more competitive 
for these problems compared to the Anderson acceleration methods. Among 
the vector extrapolation methods, the polynomial methods provide a better 
performance in general (some exceptions can be seen in the experiments 
below). 

• The main drawback of the Anderson acceleration methods concerns the nu¬ 
merical treatment of the associated minimization problem, since most of the 
difficulties come from ill-conditioning. This might be improved by including 
suitable preconditioning techniques (here the methods were implemented in 
a standard way, [7^531134] ). However, it is remarkable that when the Ander¬ 
son acceleration methods work, their performance is in general comparable 
to that of some vector extrapolation methods. 

The structure of the paper is as follows. Section [2] is devoted to a description 
of the two families of acceleration techniques considered in this study. This 
also includes some comments on the implementation and convergence results. 
The application of both techniques to the methods ([2]) is studied in Section 
E] through a plethora of numerical experiments involving the computation of 
different types of traveling waves: ground states, classical and generalized soli¬ 
tary waves as well as periodic traveling waves. The numerical study will be 
focused on the two main motivations of the paper: the improvement of the 
efficiency and the extension of application of the methods to computationally 
harder problems and where the iteration is initially not convergent. Finally, 
Section 0] completes the computational study with some illustrations of the 
application of the acceleration to the extended versions of the algorithms (12]) . 
treated in and suitable when the nonlinearity in ([1]) contains several homo¬ 
geneous terms with different degree. Some concluding remarks are in Section 

El 


2 Acceleration techniques 


Besides the local character of the convergence, in some cases and compared to 
other alternatives, fixed point algorithms has the additional disadvantage of 
a slow performance. In what follows, several techniques of acceleration will be 
considered and applied to the methods (|2]), with the aim of improving their 
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efficiency. Furthermore, as in the case of the classical algorithm, CD] , some 
cases of divergence will be transformed to convergent iterations. 

This section introduces two families of acceleration techniques: the vector 
extrapolation methods (VEM from now on) and the Anderson acceleration 
methods (AAM). We will include a description of the schemes (including some 
convergence results) and some comments on implementation. 


2.1 Vector extrapolation methods 


The hrst group of acceleration techniques consists of vector extrapolation 
methods (VEM). For a more detailed analysis and implementation of the 
methods see |25ll35ll51ll7UII41ll2UII22ll64j and references therein. Here we will de¬ 
scribe the general features of the procedures and their application to (|2]). 

Two families of VEM are typically emphasized in the literature. The hrst 
one covers the so-called polynomial methods; they include, as the most widely 
cited, the minimal polynomial extrapolation (MPE), the reduced rank extrap¬ 
olation (RRE) and the modihed minimal polynomial extrapolation (MMPE) 
methods, |25|l35|l51|l7U|l41|l66|l61j . The second family consists of the so-called 
e-algorithms; typical examples are the scalar and vector e-algorithms and the 
topological e-algorithm, [T8]l20ll41ll7T| . All the methods share of course the idea 
of introducing the extrapolation as a procedure to transform the original se¬ 
quence {un\ of the involved iterative process by some strategy. The polynomial 
methods are usually described in terms of the transformation {k < m) 


Tfc : R'" —^ R”* (5) 

^ t tjijk ^Un^k (^n,k^ Un^k^ Vnk^^n, 

where 

• A Vj^i 

• /yUn.k {i = 1, 2) denotes the m x k matrix of columns A*m„, ..., A^Un+k-i- 

• Vn,k stands for the mx k matrix of some columns ..., with V* ^ as 
the adjoint matrix of 14,fc (conjugate transpose). 

In (jSj), stands for the Moore-Penrose generalized inverse of A, dehned as 
= {A*A)~^A*, [30ll52ll38j . Different choices of the vectors jT,...,/c 
lead to the most widely used polynomial methods: 

(i) Minimal polynomial extrapolation (MPE): = Aun+j-i-,3 = 1,..., A:. 

(ii) Reduced rank extrapolation (RRE): = A^Un+j-iA = 1,..., fc. 
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(iii) Modified minimal polynomial extrapolation (MMPE): = Vj,j = 

1,... ,k, for arbitrary, fixed, linearly independent vectors Vi,.. .Vk G M™. 

The formulation of the VEM may follow an alternative approach, [6^70] . The 
transformation (j5]) can be computed in the form 


k k 

^n.k ^ ) "y 1 ) ( 6 ) 

j=0 j=0 

where the coefficients jj are obtained from the resolution (in some sense) of 
overdetermined, inconsistent systems 


fc-i 

y ) diWn+i ( 7 ) 

j=0 

for some vectors G M'”. Different methods emerge by combining dif¬ 

ferent choices of the norm where the residual vector diWn+i = Wn is 
minimized with suitable vectors Wj,Wj G M™'. Thus, for example, assuming 
k < m, we have: 

• RRE is obtained by writing ([2]) in the form 


k-l 

tn.k y ) /dj^Un-\-j 1 

j=0 

where the fdj solve ([7j) with Wj = A‘^Uj,Wj = Auj and the Euclidean norm 
with equal weights is used. 

• MPE is obtained by using (|6]) with 


7j = k' , > <0<J <k, 

Ai=o w 

where = 1 and the Cj,0 < j < k — 1 solve ([7]) with Wj = Auj,Wj = 
—Auj+k in the sense of minimization with the Euclidean norm with equal 
weights. 

• MMPE is obtained from ([6]) but where instead of ([7]) a system of the form 


fc-i 

y diQ j (^W n-\-i} Q j{Wn) 1 j 1, . . . , /c, (8) 

i=0 

is used. In (0) Qj{y) = (e^-, y) = yj, being y = {yi,..., y^)^. 

These formulations can be unihed by a representation with determinants, 
[7n|66f61|l^ . This writes the extrapolation steps tn,k in the form 
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( 9 ) 




^n+1) ■ ■ ■ 1 ^n+k) 


with 


o’o c’’! . 

^ 0,0 %,1 . UQ^k 

D((To, ..., (Tfc) — : : . : , 


( 10 ) 


^k—1,0 ^k—1,1 . ^k—l,k 

where the Uij are scalars that depend on the extrapolation method and where 
the expansion of flTU]) is in the sense 


k 

D{ao,... ,ak) = '^aiNi, ( 11 ) 

i=0 

with Ni the cofactor of cTj in the hrst row. Thus in the case of the numerator 
in ([9]), formula flTT|) is a vector, while in the case of the denominator in (p|). 
formula flTT|) is a scalar. (See [23] for the interpretation in terms of the Schur 
complement of a matrix.) The three previously mentioned polynomial methods 
correspond to the following choices of Uij: 

• AlPE. UiJ • 

• HHE. Ui^j • 

• MMPE: Uij = (ci+i, Aun+j), where ei,...efc are linearly independent vec¬ 
tors in 

A second family of VEM is called the e-algorithms. A description of them may 
start from the scalar e-algorithm of Wynn, [HITS]. This scalar extrapolation 
method can be derived from the representation (jH]), flTUD (in the scalar case) 
with Uij = Aun+j+j,* = 0,..., A; — 1; j = 0,..., A; (which are scalars in the 
scalar case). The corresponding ratio of determinants 


^n,k ('ki.^n) 


^n+l) • • • ) 'U-n+k) 


( 12 ) 


is called the classical e- (or Shanks Schmidt SS) transform, [671I68II76] . This 
ratio can be evaluated recursively for increasing k and n without the compu¬ 
tation of determinants or Schur complements. The corresponding formulation 
is 
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(13) 

(14) 


e 

e 


(n) 

-1 

(n) 

fc +1 


= 0 , 


(n) 

ey = Ur, 


= e 


(n+l) 

k-1 


+ U 


(n+l) 


n = 0,1, 2,... 


0 , 1 , 2 ,.. 




where := efc(M„), := {ek{Aun))~^ and works along diagonals on n + 
k constant. Thus, from ([12]), formulas flT3|) . flTT|) compute each entry of a 
triangular array in terms of the previous entries. 


The extension of the scalar e-algorithm to the vectorial case was carried out 
by Brezinski, [T^], and Wynn, USEZ], by using different definitions of ‘inverse’ 
of a vector, see IZD]. Wynn suggests to consider the transpose of the Moore- 
Penrose generalized inverse of a vector. 


w = (15) 

I |w| p 

leading to the vector e-algorithm (VEA), whose formulas are of the form ([T3|), 
flTTll where the scalars are substituted by vectors and flTTl) makes use of flTHll . 
This implies that the e-transform flT^ is understood in the above described 
vectorial sense. This is called the generalized Shanks Schmidt (GSS) transform, 
|18] . On the other hand, Brezinski dehnes the inverse of pair of vectors {v,w) 
such that {v,w) 7 ^ 0 as the pair of vectors where 


_i V -1 w 

w = - --, V = - -- 

\WiV) 

Thus, v~^ is called the inverse of v with respect to w and viceversa. This 
definition leads to the so-called Topological e-algorithm (TEA), when an ar¬ 
bitrary vector y is hxed and the inverses of Ae^ and are considered 

with respect to y, that is 




y 

{y, Ae^fc^+i) 


The recursive formulas are 


eLI = 0, 
4i+i = E 


(n) 

eo = Un, 


(n+l) 

2k- 


n = 0,1,2,..., 

? + (Ae£^) , k,n = 0,l,2,..., 


(n) 

4fc+2 = e 


(n+l) 

2k “T 


Ae 


(n) 

2k 


(Ae£,.A4;') 


W\ 


A:, n = 0,1, 2,... 


(16) 


Brezinski proved, [18], the connection with the GSS-transform, showing that 









4?=efcK), 4fcVi = (efcM) ^ 1 

{y,ek[un)) 

For an efficient implementation of ffT 6 |) see HU. Thns (TEA) corresponds to 
take Uij = Q{un+i+j) = {y,Un+i+j) in ([9]), 


The mechanism of working of the VEM can be described as follows, see 
[661I70II64] for details. One starts from assnming an asymptotic expression for 
the seqnence Un of the form 


Un = u + '^WjX^, n ^ oo, (17) 

i=i 

where u G and Xj G C, ordered snch that |Aj| > |Aj_|_i|,Aj 7 ^ 0,1, Aj 7 ^ 
Xj if i ^ j with only a hnite nnmber of Xj having the same modnlns. The 
expansion ffT7|l can be generalized by considering, instead of constant vectors 
Wj, polynomials Pj{n) in n with vector coefficients of the form 


Pi 


PM) = Y.^31 

1=0 



(18) 


with {vjo ,..., Vjpj} linearly independent in MP and where if |Aj| = |Aj+i| then 
Pj > Pj+ij [SS]- For simplicity, the description below will make nse of (fT7|) . 


The asymptotic expansion 0171) is considered in a general vector space (hnite 
or inhnite dimensional) where the iteration is dehned. It is nnderstood in the 
sense that any trnncation differs from Un in less than some power of the next 
A. This means that for any positive integer N there are K > 0 and a positive 
integer uq that only depend on N snch that for every n > hq 


N-l 

\\Un-U- 

In particnlar, the case N = 1 allows to identify u as limit or anti-limit of the 
seqnence according to the size of Ai. (That is, if |Ai| < 1 then lim^^ooWn 
exists and eqnals u. If |Ai| > 1 then hm„_j.ooM„ does not exist and u is called 
the anti-limit of the seqnence m„.) Under these conditions, several resnlts of 
convergence for MPE, RRE, MMPE and TEA are obtained in the literatnre, 
[IFIllTn] and references therein. For these methods, one can End an extrapolation 
step K snch that 

I \tn,K -u\\ = 0(A"+i), 00. (19) 
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The estimate ffT^ may explain the convergent behaviour of the extrapolation 
in some cases. If the A’s are identified as the eigenvalues of the linearization 
operator of the iteration at the limit (or anti-limit) u, then the extrapolation 
has the effect of translating the behaviour of the iteration to an eigenvalue 
Ak +1 that may be into the unit disk, even if the previous ones are out of it. 
Hence, u may be anti-limit for the original iteration and the extrapolation 
converges to it. These results are extended to the defective linear case with 
more general polynomials flT8|l in |65] . 

The integer k is related to the concept of minimal polynomial T’(A) of a matrix 
A with respect to a vector v, |39]I41II70] : this is the unique polynomial of least 
degree such that 

P{A)v = 0. 

Thus in the case of linear iteration with matrix A, k is taken to be the degree 
of the minimal polynomial of A with respect to the hrst iteration uq. In the 
case of a nonlinear system written in hxed point form 


a: = P{x), (20) 

then K is theoretically defined as the degree of the minimal polynomial of 
A = P'{u*) with respect to Mq, where u* is a solution of (120]) . In contrast with 
the linear case, there is no way to determine k in advance for the nonlinear case. 
This forces to consider several strategies for the choice and the corresponding 
implementation, see the discussion in cni and the comments here below. 

We also mention that in the linear case, the extrapolation methods MPE 
and RRE are mathematically equivalent to the method of Arnold!, [00], and 
the GMRES, [00], respectively, see [02], while the MMPE is mathematically 
equivalent to the Hessenberg method, [00] and TEA to the method of Lanczos, 
[^, see [02IITT] . 

Efficient and stable implementation of the RRE, MPE and MMPE methods 
by using QR and LU factorizations can be seen in [OOllTO] . For the case of 
TEA, see im and [2TII0T] for VEA. The implementation is usually carried out 
in a cycling mode. A cycle of the iteration is performed by the following steps 
: consider a method ([2]) with a stabilizing factor s satisfying (PI), (P2). Given 
Mo 7^ 0 and a width of extrapolation mw > 1, for I = 0,1,..., the advance 
I i —y / -)-1 is! 

(A) Set to = ui and compute mw steps of the fixed-point algorithm: 

Ltn+i = s{tn)N{tn),n = 0 ,... mw - 1 . 

(B) Gompute the extrapolation steps ([9]) with any of the methods described 
above and n = 0, ..., mw. 
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(C) Set Ui+i = tmw,i,to = Ui+i and go to step (A). 


The cycle (A)-(B)-(C) is repeated until the error (residual or between two 
consecutive iterations) is below a prehxed tolerance, a maximum number of 
iterations is attained or the discrepancy between the stabilizing factor at the 
iterations and one is below a prehxed tolerance. 


The width of extrapolation mw depends on the choice of the technique: mw = 
K +1 for MPE, RRE or MMPE and mw = 2k for VEA or TEA, [70]. Since k is 
generally unknown, some strategy for the implementation must be adopted. In 
practice, as discussed in ra, the methods are implemented with some (small) 
values of mw and take that with the best performance. It may be also different 
for each cycle, although quadratic convergence is not expected if k is too small. 
This choice of k will be computationally studied in some examples in Section 
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The hypotheses for the expansion fITTD include the conditions Xj 7 ^ 1 Vj. In 
many problems for traveling wave generation, A = 1 appears as eigenvalue 
of the iterative technique (of hxed point type) although under especial cir¬ 
cumstances that allow to extend the convergence results in some sense. This 
especial situation is related to the presence of symmetries in the equations 
for traveling waves. In order to extend the results of convergence to this case, 
one has to consider the orbits by the symmetry group of the equations and 
interpret the convergence in the orbital sense, that is a convergence for the 
orbits. The description of this orbital convergence can be seen in [3]. 


Finally, local quadratic convergence is proved in ca (see also 
for the four methods and VEA for a general nonlinear system fl2U]) and under 
the following hypotheses on J^: 


• The Jacobian matrix does not have A = 1 as eigenvalue. 

• k is the degree of the minimal polynomial of with respect to Uq — u*. 

• the algorithm is implemented in the cycling mode where k is chosen in the 
Ath cycle as the degree of the minimal polynomial of J^'{u*) with respect 
to ti-\ — u*. 


This result can be extended to the case where admits a i/-parameter {v > 1) 
group of symmetries and, consequently, A = 1 is an eigenvalue of J^'(n*), 
by using the reduced system for the orbits of the group and in the orbital 
sense, [3]. This would lead to quadratic orbital convergence but linear local 
convergence. 
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2.2 Anderson acceleration methods 


A second family of acceleration techniques considered here is the so-called 
Anderson family or Anderson mixing, |6] . It is widely used in electronic struc¬ 
ture computations and only recently it has been analyzed in a more general 
context, [7311531182] ■ (To our knowledge, this is the hrst time that AAM are 
applied to accelerate traveling wave computations.) The main goal of the ap¬ 
proach consists of combining the iteration with a minimization problem for 
the residual at each step. For linear problems, this technique is essentially 
equivalent to the GMRES method, [60)1301138] . The stages for an iteration step 
are as follows: Given uq ^ 0, nw > 1, set Lu\ = s{uo)N{uo). For k = 1,2,... 


• Set Uk = min{nt<;, k} 

• Set Fk = (/fc-Hfc, • • •, /fc) where fi = Lui - s{ui)N{ui) 

• Determine = (aQ^\ ..., that solves 


rik 

min \\Fka\\, ^^* = 1 


• Set 


Luk+i — 


Y,(^i^'’s{uk-nk+i)N{uk -Uk+i) 
i=0 


( 21 ) 


The resolution of the optimization problem fl2ip is the source of the additional 
computational work of the acceleration. One way to reduce this extra effort is 
the so-called multisecant updating [MIISG] . (This also clarihes the connection 
with quasi-Newton methods.) This technique consists of writing the problem 
in an equivalent form 


min 11/- J'fcyll, J'fc = (A/fc_„^,..., A/fc_i), A/i =/i+i -/i, (22) 

7=(70r-->7nfc) 

but with a more direct resolution, and determining the acceleration from it. 
The general step becomes: 

• Set Hk = min{ntc, k} 

• Determine 7 *^^^ = ( 7 o*^\ ..., 7 );^^^) by solving fl22]) . 

• Set 


(k) 

«o 


(fc) 

= 7o b 


a. 


(k) 


(k) 

= 11 -1: 


{k) 


I <i <nk-l, = 1 - 


(fc) 

Ink-l 


• Set 


Luk+i — 


Y,(^\’'^s{Uk-nk+i)N{Uk -Uk+i)- 
i=0 
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As mentioned in m, if J^k is full-rank, the solution of the minimization prob¬ 
lem can be written as 7 *^^^ = Fk)~^fk and the Anderson acceleration 

has the alternative form 


LUk+l LUk Gkfki 

Gk = —I+ ('Hfc -I- Fk){J^k ^k)~^^k 1 'Hk = {Auk-nik, ■ ■ ■ y ^Uk-i), 

tli. 

(Note that Gk can be viewed as an approximate inverse of the Jacobian of 
Lx—N{x)). The formulation fl25D motivates the generalization of the Anderson 
mixing, [31]. This is performed replacing (J^J^k)~^J^k in 11231) by some Vk € 
l^nxm satisfying 



= /, 


and fl23]l becomes 


LUk+l LXk Gkfki 

Gfc = —/-1- ("Hfc +Hk = i^Uk-rriki ■ ■ ■ 1 ^Uk-l). (24) 

The resulting methods are collected in the so-called Anderson’s family, [751131] . 
Two particular members are emphasized: the Type-I method (denoted by AA- 
I from now on), which corresponds to Vk = {'HfJ^k)^^'Hk in fl2l|) and Type-II 
method (or AA-II), which is the original Anderson mixing 1I23]1 . 

To our knowledge, some convergence results can be seen in [7311571172] . In [73] 
the authors identify some Anderson methods for linear problems and in some 
sense with the GMRES method and the Arnoldi (FOM) method; some conver¬ 
gence results can be derived from this identification. In [37] the equivalence 
with GMRES for linear problems is completely characterized. Finally, [72] 
gives some proofs of convergence of the Anderson acceleration when applied 
to contractive mappings: g-linear convergence of the residual for linear prob¬ 
lems under certain conditions when nw = 1 and local r-linear convergence in 
the nonlinear case. (These types of convergence are defined in the paper.) 

On the other hand, as observed in ca, the implementation of AAM should 
be carried out by attending to three main points: a convenient formulation 
of the minimization problem fl?T]) . a numerical method for its efficient resolu¬ 
tion and, hnally, the parameter nw, which plays a similar role to that of the 
extrapolation width mw in the VEM. In our computations below, we have 
followed the treatment described in m- This is based on the use of the un¬ 
constrained form (122|) and its numerical resolution with QR decomposition. 
For other alternatives in both problems, see the discussion in [731134] and the 
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references cited there. (According to our results below, the use of alternative 
preconditioning techniques might be recommendable in some cases.) As far as 
the choice of nw is concerned, a similar strategy to that of mw will be used, 
since our experiments, na, suggest that nw (as mw) strongly depends on the 
problem under study and large values are not recommended. Finally the codes 
are implemented by retaining the dehnition of Uk = min{nta, k} since other 
alternatives, [82], did not improve the results in a relevant way. 


3 Numerical comparisons 


Presented here is a comparative study on the use of VEM and AAM as accel¬ 
eration techniques from the Petviashvili type methods ([2]) in traveling wave 
computations. The comparison is organized according to two main points: the 
type of traveling wave to be generated and the elements of each family of 
methods to be used for the generation. The hrst point includes the following 
case studies: 


( 1 ) 

( 2 ) 

(3) 


Classical solitary waves, generalized solitary waves and periodic traveling 
waves of the four-parameter Boussinesq system, [M. 

Localized ground state solutions of NLS type equations, iHHoi. 

Highly oscillatory solitary waves of the one- and two-dimensional Ben¬ 
jamin equation. 


This plethora of waves attempts to discuss and overcome different computa¬ 
tional difficulties and with the aim of establishing as more general conclusions 
as possible. As for the second point, each family of techniques has been rep¬ 
resented by the following methods: 


• MPE and RRE standing for polynomial extrapolation methods. 

• VEA and TEA standing for e algorithms. 

• AA-1 and AA-11 standing for the AAM. 

For simplicity, the acceleration will be applied to the Petviashvili method l|2|), 
l]T|) with 7 = p/{p — 1). Due to the similar behaviour of the methods of the 
family ([2]), illustrated in [3], the conclusions from the corresponding results 
can reasonable serve when the Petviashvili method is substituted by any of 

O- 


In all the cases considered, the traveling wave prohles appeared as solutions 
of initial value problems of ode’s. Their discretization to generate approx¬ 
imations to the prohles, was carried out in a common and standard way. 
The corresponding initial periodic boundary value problem (on a sufficiently 
long interval) was discretized by using Fourier collocation techniques, [TTII^ : 
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the discretization leads to a nonlinear system of algebraic equations for the 
approximate values of the profile at the collocation points or for the discrete 
Fourier coefficients of the approximation. This system is iteratively solved with 
the classical Petviashvili method ([2]), (jl]) along with the selected acceleration 
technique. This will be described in each equation considered below. 

Several stopping criteria for the iterations are implemented: 

• A maximum number of iterations. 

• The iteration stops when one of the following quantities are below a prefixed, 
small tolerance TOL: 

(i) The difference in Euclidean norm between two consecutive iterations 

I l^n+1 ||) ^ 0,1,... 

(ii) The residual error (also in Euclidean norm) 

RESn = \\Lun - N{un)\\, n = 0, 1 ,... ( 25 ) 

(iii) The discrepancy between the stabilizing factor and (in case of conver¬ 
gence) its limit one 

SFEn = \s{un) - l\, n = 0,l,... ( 26 ) 

The numerical experiments that form the comparative study are of different 
type: 

• Eor several values of k (in the case of VEM) and nw (in the case of AAM) 
we have computed the number of iterations required by each method to 
achieve a residual error below the prehxed tolerance. This allows to compare 
some performance of the methods between the two families, among different 
techniques within a same family and indeed with the Petviashvili method 
without acceleration. 

• Some eigenvalues of the iteration matrices for the classical hxed point al¬ 
gorithm and the Petviashvili method have been computed (with the corre¬ 
sponding standard MATLAB function) in order to explain the behaviour of 
the second one, [3] and how the acceleration eventually changes it. 

• The form of the approximate prohles and some experiments to check their 
accuracy are also displayed. 


3.1 Traveling wave solutions of Boussinesq systems 


In this first example we study the numerical generation of traveling wave 
solutions of the four-parameter family of Boussinesq system 
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( 27 ) 


7]t "1“ Vjx “1“ “t“ OjUxxx bvjxxt 0, 

Uf -\- Tjx UVjx CTjxxx d^xxt 0 ; 

where rj = //(x,t),ti = u{x^t),x G M,t > 0 and the four parameters a.b^c^d 
satisfy 


o + 6 = i(«"-i). c + <i = i(l-e^), (28) 

with some constant 9'^ e [0,1], [TB]llllll2j . System fE7|) appears as one of the 
alternatives to model the bidirectional propagation of the irrotational free snr- 
face flow of an incompressible, inviscid flnid in a nniform horizontal channel 
nnder the effects of gravity when the snrface tension and cross-channel varia¬ 
tions of the flnid are assnmed to be negligible. If ho denotes the nndisturbed 
water depth, then r]{x,t) stands for the deviation of the free snrface from hg 
at the point x and time f, while u{x,t) is the horizontal velocity of the flnid 
at X and at the height y = dh^ (where j/ = 0 corresponds to the channel 
bottom) at time t. For the derivation of (12711 from the two-dimensional Euler 
equations and the mathematical theory see mug. For the modihcation of 
(12511 to include the influence of surface tension see [2511271128) . 

The Boussinesq system (127)1 admits different types of traveling wave solutions. 
First, being an approximation to the corresponding two-dimensional Euler 
equations in the theory of surface waves, it is expected to have classical solitary 
wave solutions. They are solutions of the initial value problem of (127)1 of smooth 
traveling wave form rj = r]{x — Cst),u = u{x — Cgt) with some speed > 0 and 
decaying to zero as X = x — Cgt ^ ±cxo. Substitution into (127)1 and after one 
integration the prohles u = u{X),r] = rj^X) must satisfy the ode system 


Cs(l — hdxx) —(1 + h 

(1 -I- cdxx) Cs(l — ddxx) j yw 

The problems of existence, asymptotic decay and stability of solutions of (12911 
have been analyzed in many references and for particular values of a, b, c, d, 
see [55] and references therein. Furthermore, in the same reference, linearly 
well-posed systems (129)1 may be studied as a hrst order ode system and, 
based on normal form theory, a discussion on the values of the parameters 
leading to Boussinesq systems admitting solitary wave solutions with speed 
Cs > 1 is established. According to it, two classes of systems can be distin¬ 
guished. The hrst one admits classical (in the sense above dehned) solitary 
wave solutions. This group contains the Bona-Smith system {a = 0,b = d = 
(39^ - l)/6,c = (2 - 302)/3,2/3 < 9^ < 1), [I5], or the BBM-BBM system 
(a = c = 0,b = d= 1/6), [10]. (The classical Boussinesq system, which corre¬ 
sponds to a = 6 = c = 0, d = 1/3, is also in this group, although it is out of 


ur] 
\ 2 


(29) 
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the general discussion of [33] and has been studied separately, [HI].) A second 
class of Boussinesq systems admits generalized solitary wave solutions, that is 
traveling wave prohles which are not homoclinic to zero at inhnity but to small 
amplitude periodic waves, [HD]- The KdV-KdV system {a = c> 0,b = d = 0), 
[I3II3, is an example of this second group. 

Finally, the existence of periodic traveling wave solutions of fl27|) is studied in 
[26] by applying topological degree theory for positive operator to the corre¬ 
sponding periodic initial value problem posed on an interval {—1,1) and some 
cnoidal wave solutions of the BBM-BBM system are computed. A smooth 
periodic traveling wave solution rj = ri{x.Cst),u = u^x.Cgt) with some speed 
Cs > 0 must satisfy the ode system 


Cs(l bdxx) (It ^dxx) 
y (It cdxx) Cs(l ddxx) 



^ur]^ 





[K 2 ) 


(30) 


for some real constants Ki, K 2 . These are related to the period parameter 1. 
The resolution of fl30|l involves a modihed system for which these constants 
of integration are set to zero, [26]. This is briefly described as follows. One 
hrst searches for constant solutions rj = Ci,u = C 2 of fl30|l . This leads to the 
system 


/ 

V 


-1 




^C^C2\ 

Uj 


\K2) 


5 


which can be solved as a cubic equation for C 2 '- 


Y - T (c2 - l-f- Ki)C2 -K 2 - CsK, = 0, (31) 

Oi = c.Oa - Ci/2 - K 2 . (32) 

Once ( 01 , 02 ) is obtained (there may be more than one solution indeed) the 
differences rj = t] — Ci,u = u — C 2 must satisfy 


^Cs{l - bdxx) - O2 
~(1 T cdxx) 


— (1 -|- adxx) — Cl 

c^(l - ddxx) - O 2 




f~~\ 


V 

urj 1 



72 

1 

[uj 

Vt/ 


(33) 


This strategy will be considered in the numerical generation of the prohles 
in (l30|) : the system fl3^ wil be discretized to compute approximations to the 
variables rj, u and to the variables rj = rj + Ci,u = u + C 2 from them. 
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In order to generate numerically classical and generalized solitary wave solu¬ 
tions of (HT} the corresponding periodic value problem of fl2^ on a long enough 
interval {—1,1) is discretized with a Fourier collocation method leading to a 
discrete system of the form 


+ Cs{In.-dD^) 

- -v' 

L 



^Uh-Vh 

Uh-' 

\ 2 


N{r)h,Uh) 


(34) 


where rjh.Uh G are approximations T]h,j ~ ~ u{xj) to the values 

of a solution of fl30|l at the grid points Xj = —l+jh, h = 2l/m,j = 0,.. .m —1, 
D is the pseudospectral differentiation matrix, is the mxm identity 

matrix and the nonlinear term N, which is homogeneous of degree p = 2, 
involves Hadamard products. In the case of periodic traveling waves and as 
was mentioned above, system fIMD will be substituted in the implementation 
by 


^Cs{Im - bD^) - C^Im -{Im + aD^) - CJm\ f dh 

^ + Cs{Im - dD‘^) - C2lm) \uh 

for the approximations rjh, Uh to the fj, u variables at the grid points and where 
(71,(72 are previously known from the resolution of fl3TH . (13^ . Then rjh = 
Vh + Cl, Uh = Uh + C2- 

The methods ([2]) along with the corresponding acceleration technique are then 
applied to the discrete systems (l3T|) and ([35]). The implementation is performed 
in the Fourier space; for example (IMD becomes 


Uh-Vh 


(35) 


7(i + 6(y)5 -(i-o(y)7 



Thus the 2m x 2m system (l3T|) is divided into m blocks of 2 x 2 systems for 
the corresponding p-th discrete Fourier coefficients , {uh)p , —m/2 < p < 
m/2. (For simplicity, we assume that m = 2^ for some s > 1.) Alternatively, 
fl3Tll can be written in the form 
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= Th 

'’"‘1 = 1 

V“V 




Ah * {uh.r]h) + Bh* 
Bh * {uh-r]h) + Ch* 



where * denotes periodic convolution and if a; = exp(—27rf/m), the vectors 
Ah, Bh, Ch have discrete Fourier coefficients 


{Ah)p — 
{Ch)p = 


1 - a iff 

A(p) ’ 
A(p) ’ 


{Bh)p — 
{Dh)p = 


c.(i+b{^y) 

A(p) 

c.(l + d(lA)^) 

A{p) 


A(p) = + )(i + ^ 


m m 

— < p < —. 
2 “ 2 



In order to explain the behaviour of the iteration, the size of the eigenvalues 
of the iteration matrix will be relevant in the numerical study. In this case, 
the corresponding iteration matrix of the classical fixed point iteration at a 
solution u* = {7]y,uD has the form 


^ ^ diag(M^) diag(r/;^)\ 

V 0 diag(M^) ) 

(where diag(n) stands for the diagonal matrix with diagonal entries given 
by the components of n G M™'). Some information on the spectrum of S is 
known. We already have the eigenvalue A = 2, corresponding to the degree of 
homogeneity of the nonlinear part, with u* = {Ph^'^h) eigenvector. Also, 

the application of D to fIMl) leads to 


CsDpl - Dul = D{plul) = ul.Dpl + pI-DuI 
-Dill + c.(l- Dnl = D = ulDul, 


which means that A = 1 is an eigenvalue of S and {Dph, B>Uh)^ is an associated 
eigenvector. This corresponds to the ‘translational’ invariance of 0271) . 


Three particular systems of 0271) will be taken to illustrate the numerical gen¬ 
eration of traveling waves. The first one is the classical Boussinesq system 
{a = b = c = 0,d = 1/3), 
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( 36 ) 


rit + Ux + {w)x = 0 , 

1 

“1“ Vx UUx '^Uxxt O' 

which is known to have classical solitary wave solutions, [M]' The second one 
is the so-called KdV-KdV system (a = c= l/6,b = d = 0) 


dt T Ux -|- {^u)x T p^xxx 0, 

6 

T ^x T 'uUx ~'f]xxx 0, (07) 

6 

that admits generalized solitary wave solutions, [mi]. Finally, in order to 
illustrate the numerical generation of periodic traveling waves, the BBM-BBM 
system {a = c = 0,h = d = 1 / 6 ), 


Vt + Ux + {vu)x - -Vxxt = 0, 

6 

1 , , 

Ut "F Tjx -|- UUx ~^^xxt 0 , (08) 

6 

will be taken, [26]. 


3.1.1 Numerical generation of classical solitary waves of 

In the case of system fl36l) a hrst experiment of comparison of the acceleration 
techniques has been made by taking c^, = 1.3 and a hyperbolic secant prohle as 
initial iteration with I = 64 and m = 1024 collocation points. The Petviashvili 
method (Ej), ( 0 ]) with 7 = 2 was hrst run, generating approximate rj and u pro- 
hles as shown in Figures [I](a) and (b) while Figures [T](c) and (d) stand for the 
corresponding phase portraits of the approximate prohles in (a) and (b). (They 
show the classical character of the solitary waves, represented as homoclinic 
to zero orbits with exponential decay, [54]. ) The accuracy of the iteration is 
checked in Figure O Figure [2](a) illustrates the convergence of the sequence 
Sn = s{7]n, Un) of Stabilizing factors, computed with the corresponding formula 
dl]) and the optimal value 7 = 2 . The discrepancy fl26]l is below the tolerance 
TOL = 10“^^ in 77 , = 62 iterations, while the hrst residual error below TOL is 
9.09248977 — 14 at n = 76. (This also happens in the rest of the experiments: 
when the procedure is convergent, the error |1 — s„| achieves the tolerance 
before the residual error; therefore, the control on this last one is a harder test 
and will be adopted as the main one to stop the iteration.) Convergence is 
also conhrmed by Table [H This shows the six largest magnitude eigenvalues 
of the iteration matrix (l36]) (hrst column) and of the iteration matrix (at the 
same iterate) of the Petviashvili procedure (second column) both at the last 
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(a) 


(b) 



T1 U 


(c) 


(d) 


Fig. 1. Approximate profiles generated by the Petviashvili method ([2]), ([1]) for (j36l] 
with Cs = 1.3: (a) ij, (b) tt; (c) Phase portrait of r], (d) Phase portrait of u. 



(a) 


(b) 


Fig. 2. Convergence results of the Petviashvili method for (|36p : (a) Discrepancy (1261) 
for the stabilizing factor Sn = s{T]n,Un) vs number of iterations: (b) Residual error 
(j25p vs number of iterations. (Semi-logarithm scale in both cases.) 


21 













Iteration matrix S{r]f,Uf) 

Iteration matrix F'{rif,Uf) 

1.999999E + 00 

9.999999A - 01 

9.999999A - 01 

6.763242A - 01 

6.763242A - 01 

5.411229A - 01 

5.411229A - 01 

4.820667A - 01 

4.820667A - 01 

4.567337A - 01 

4.567337A - 01 

4.465122A - 01 


Table 1 

Classical solitary wave generation of (|36|) . Six largest magnitude eigenvalues of the 
approximated iteration matrix S = L~^N'[rjj, uj) (first column) and of the iteration 
matrix F'{r]f,Uf), generated by the Petviashvili method ([2]), ([3]) with 7 = 2, both 
evaluated at the last computed iterate {rjf,Uf). 

computed iterate {T]f,Uf). The first column reveals the dominant eigenvalues 
Ai = 2, A 2 = 1, both simple, while the rest is below one. The filtering effect of 
the Petviashvili method, |3], is observed in the second column; the dominant 
eigenvalue is filtered to zero (recall that 7 = 2) and the rest is preserved. Since 
A 2 = 1 corresponds to the translational symmetry of ([27j), this guarantees the 
local convergence of the method (also in the orbital sense mentioned above). 

The improvement of the performance of the Petviashvili method with several 
acceleration techniques is now computationally analyzed. A first point to study 
is the choice of the parameters k (for the VEM) and nw (for the AAM). Table 
|2]shows, for values of k between one and ten, the number of iterations required 
by MPE, RRE, VEA and TEA to achieve a residual error below TOL = 10“^^. 
(The residual error, corresponding to the last iteration is in parenthesis for 
each computation.) From these results, the following comments can be made: 

(a) For K > 2, all the methods improve the performance of the Petviashvili 
method without acceleration (cf. Figure [2](b)). The reduction in the num¬ 
ber of iterations varies in a range 50 — 70%. 

(b) In general, polynomial methods (MPE and RRE, which essentially be¬ 
haves in an equivalent way) are more efficient than e-algorithms (with 
VEA slightly better than TEA). In the best cases, the improvement is 
about 70% in the case of MPE and RRE and RRE, about 65% with 
respect to VEA and about 54% in the case of TEA. (However, one has 
to take into account that the cycle in the case of polynomial methods is 
mw = K -|- 1 and in the case of e-algorithms is mw = 2k; cf. Figure El) 

In the case of the AAM, the corresponding results are in Table El Now, the role 
of the parameter k (or mw) is played by nw. The results show that the per¬ 
formance of the methods is essentially the same. The best results are obtained 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

1 

269 

99 

631 

408 


(8.9136E - 14) 

(9.1312E - 14) 

(9.7920E - 14) 

(7.3356E - 14) 

2 

64 

48 

43 

43 


(7.4794E - 14) 

(8.3903E - 14) 

(7.7841E - 14) 

(8.1979E - 14) 

3 

43 

43 

38 

42 


(7.5001E - 14) 

(8.0682E - 14) 

(7.7395E - 14) 

(9.6255E - 14) 

4 

33 

33 

33 

37 


(7.9601E - 14) 

(8.2823E - 14) 

(7.7824E - 14) 

(8.0527E - 14) 

5 

28 

28 

31 

39 


(8.5285E - 14) 

(9.8557E - 14) 

(8.7444E - 14) 

(9.3163E; - 14) 

6 

26 

26 

29 

35 


(8.3589E - 14) 

(7.7189E - 14) 

(7.3281E - 14) 

(7.9462E; - 14) 

7 

27 

27 

33 

35 


(8.6403E - 14) 

(8.1151E - 14) 

(7.2215E - 14) 

(8.4479E - 14) 

8 

27 

25 

29 

37 


(9.7379E - 14) 

(8.0842E - 14) 

(7.3955E - 14) 

(7.4142E - 14) 

9 

23 

24 

30 

41 


(9.5276E - 14) 

(8.3798E - 14) 

(8.7013E - 14) 

(7.2068E - 14) 

10 

25 

25 

27 

35 


(7.6433E - 14) 

(8.3980E - 14) 

(9.3658E - 14) 

(8.5795E - 14) 


Table 2 

Classical solitary wave generation of (I36p . Number of iterations required by MPE, 
RRE, VEA and TEA as function of k to achieve a residual error below TOL = 
10“^^. The residual error at the last computed iterate is in parenthesis. Without 
acceleration, the Petviashvili method ([2|), Q with 7 = 2 requires n = 76 iterations 
with a residual error 9.092577 — 14. 

with nw = 8 in the case of AA-I and nw = 9 for the AA-II. On the other 
hand, as mentioned in |53ll75] , the value of nw cannot be too large, because of 
ill-conditioning. In this example, this was observed for AA-I when nw = 9 , 10. 
(The corresponding results in Table [3] were obtained by using standard precon¬ 
ditioning.) Finally, compared to the Petviashvili method without acceleration, 
the reduction in the number of iterations is in range of 50 — 80%. 

Since the implementation of the methods is different, a comparison between 
VEM and AAM should take into account several efficiency indicators. In our 
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nw 

AA-I(ntc) 

AA-II(nTi;) 

1 

38(8.4014A - 14) 

35(4.8504A - 14) 

2 

28(4.9835A - 14) 

26(5.6978A - 14) 

3 

28(5.5678A - 14) 

25(6.2897A - 14) 

4 

27(3.8773A - 14) 

22(1.4624A - 14) 

5 

22(4.4004A - 14) 

20(6.5530A - 14) 

6 

21(7.9925A - 14) 

21(2.4615A - 14) 

7 

21(2.3111A - 14) 

20(5.3227A - 14) 

8 

20(8.0666A - 14) 

20(2.7701A - 14) 

9 

20(4.5873A - 14) 

19(9.6556A - 14) 

10 

20(2.7208A - 14) 

19(6.8255A - 14) 


Table 3 

Classical solitary wave generation of (|36p . Number of iterations required by AA-I 
and AA-II as function of nw to achieve a residual error below TOL = 10“^^. The 
residual error at the last computed iterate is in parenthesis. Without acceleration, 
the Petviashvili method ([2]), Q with 7 = 2 requires n = 76 iterations with a residual 
error 9.0925E — 14. 

example, we have measured the performance by computing the residual error 
as function of the number of iterations (i. e. comparing the best results of 
Tables [2] and [3]) and as function of the computational time. The comparison 
of the methods in terms of the number of iterations is illustrated in Figure 
[3]^a). This shows, in semilogarithmic scale, the residual error as function of the 
number of iterations for the Petviashvili method without acceleration (solid 
line) and accelerated with the six selected techniques, implemented with the 
values of k and nw that, according to Tables [2] and [3], lead to the best number of 
iterations. For this example, the AA-1(8) and AA-11(9) give, for a tolerance of 
TOL = 10“^^ in the residual error, a slightly smaller number of iterations than 
the (mostly equivalent) RRE(9) and MPE(9). The initially worse performance 
of VEA(IO) and TEA(6) is corrected after the hrst cycle. Eor example, in the 
case of VEA(IO), after this hrst cycle. Figure [3](a) shows that the reduction 
in the residual error is the fastest. 

A second comment concerns the computational efficiency. Figure [3](b) shows 
(again in semi-log scale) the residual error as function of the CPU time in 
seconds for the four VEM. According to this, VEA(IO) is the most efficient, 
followed by MPE(9), TEA(6) and RRE(9). The comparison in CPU time 
of this last one (the worst one among the VEM in efficiency) with AA-1(8) 
and AA-11(9) is shown in Figure [3](c) and reveals the poor performance in 
computational time as the main drawback of the AAM for this case, (see the 
formulation and implementation described in Section |2] to attempt to give an 
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(b) 



Fig. 3. Convergence results of the Petviashvili method for (I36|) : (a) Residual error 
(j25p as function of number of iterations for the Petviashvili method without ac¬ 
celeration and for six acceleration techiS^ues with the best parameters k and nw 
(according to Tables [5] and [6]) . (b) Residual error (|2^ as function of CPU time (in 
seconds) for six acceleration techniques with the best parameters k and nw (ac- 


























Iteration matrix S{gf,Uf) 

Iteration matrix F'{gf,Uf) 

1.999999.F + 00 

1.000000.F + 00 

l.OOOOOO.F + 00 

5.625613.F - 01 

5.625613.F - 01 

-3.525656F; - 01 

-3.525656F - 01 

-3.521308F; - 01 

-3.521308F; - 01 

3.069304.F - 01 + i6.906434F; - 02 

3.069304.F - 01 + i6.906434F - 02 

3.069304F; - 01 - ii6.906434F; - 02 


Table 4 

Generalized solitary wave generation of (l37|) . Six largest magnitude eigenvalues 
of the approximated iteration matrix S = L~^N'{r]f,Uf ) (first column) and of 
the iteration matrix F'{rif,Uf), generated by the Petviashvili method ([2]), Q with 
7 = 2, both evaluated at the last computed iterate {r]f,Uf). 


explanation of it. 


3.1.2 Numerical generation of generalized solitary waves of 

Here we show the resnlts concerning the generation of approximate generalized 
solitary waves of the KdV-KdV system (|37ll . In this case we have considered 
a speed Cg = 1.3 and a Ganssian-type prohle as initial gness for rj and u, 
with I = 64 and m = 1024 Fonrier collocation points. The approximate r] 
and u prohles generated by the Petviashvili method (without acceleration) 
are displayed in Figures |l](a) and (b) respectively (observe the oscillatory 
ripples to the left and right of the main pulse), and the performance of the 
method (measured in terms of the convergence of the stabilizing factor and the 
behaviour of the residual error as function of number of iterations) is shown 
in Figures 111(c) and (d) respectively. The method achieves a residual error of 
1.150546F1 — 12 in n = 47 iterations and 7.85942277 — 14 in n = 52 iterations. 
In this case, the corresponding phase portraits in Figures [5](a) and (b) show 
the generalized character of the waves, with orbits that are homoclinic to 
small amplitude periodic oscillations at inhnity. The last computed iterate, 
corresponding to this residual error, is used to evaluate the iteration matrices 
S{rif,Uf) and F'{r]f,Uf) of the classical hxed point and Petviashvili method, 
respectively. The associated six largest magnitude eigenvalues are shown in 
Table H (The generalized character of the computed solitary wave is also 
noticed by the presence of conjugate complex eigenvalues in the linearization 
matrix at the wave, cf. Tabled]) 

The performance of the acceleration techniques is hrst checked in Tables [5] and 
|6] (respectively) which are the analogous to Tables [2] and [3] respectively for the 
generalized solitary wave generation. The conclusions are the same as those of 
the generation of approximate classical solitary wave prohles for system fl36ll : 
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(a) (b) 




(c) 


(d) 


Fig. 4. Generalized solitary wave generation of (I37h . (a)-(b) Approximate rj and u 
profiles generated by the Petviashvili method ([2]) , (H]) for ([37)1 ; (c) Discrepancy (l26l) 
for the stabilizing factor Sn = s{r]n,Un) vs number of iterations; (d) Residual error 
()25ll vs number of iterations. (Semi-logarithm scale in both cases.) 




(a) (b) 

Fig. 5. Generalized solitary wave generation of (1371) . (a) Phase portrait of rj, (b) 
Phase portrait of u. 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

1 

93 

64 

283 

67 


(8.2739E - 14) 

(9.4330E - 14) 

(6.7393E - 14) 

(6.5919E - 14) 

2 

42 

43 

81 

49 


(9.3934E - 14) 

(8.6837E - 14) 

(6.1237E - 14) 

(8.7121E - 14) 

3 

37 

37 

38 

37 


(9.8740E - 14) 

(5.6549E - 14) 

(6.7020E - 14) 

(7.5320E - 14) 

4 

33 

32 

39 

30 


(7.5531E - 14) 

(6.9363E - 14) 

(3.6501E - 14) 

(4.2863E - 14) 

5 

28 

28 

32 

30 


(6.6301E - 14) 

(7.0496E - 14) 

(5.9777E - 14) 

(8.5244E; - 14) 

6 

28 

28 

29 

28 


(7.8419E - 14) 

(7.4939E - 14) 

(3.2861E - 14) 

{7.7335E - 14) 

7 

27 

27 

32 

32 


(3.0803E - 14) 

(3.2406E - 14) 

(6.5803E - 14) 

(4.9439E - 14) 

8 

23 

23 

28 

29 


(9.3872E - 14) 

(8.8601E - 14) 

(8.0707E - 14) 

(9.1792E - 14) 

9 

23 

23 

26 

27 


(3.0380E - 14) 

(3.1758E - 14) 

(6.4506E - 14) 

(8.0379E - 14) 

10 

24 

24 

24 

28 


(5.8031E - 14) 

(5.0834E - 14) 

(7.1291E - 14) 

(5.8663E - 14) 


Table 5 

Generalized solitary wave generation of (j37p . Number of iterations required by 
MPE, RRE, VEA and TEA as function of k to achieve a residual error below 
TOL = 10“^^. The residual error at the last computed iterate is in parenthesis. 
Without acceleration, the Petviashvili method (l2|), (jH) with 7 = 2 requires n = 52 
iterations with a residual error 7.8594E — 14. 


in terms of the number of iterations, AAM give the best performance and 
amongst the VEM, the extrapolation methods MPE and RRE are (in this case 
slightly) more efficient than the vector e-algorithms VEA and TEA, see Figure 
[6](a). The ranking is the opposite when residual error is measured in terms 
of the computational time. Figure |6](b) shows that VEA is the fastest and 
MPE the slowest. Even though, this is much faster than any of the Anderson 
algorithms, as shown in Figure [6](c). 












(b) 



Fig. 6. Convergence results of the Petviashvili method for (1371) : (a) Residual error 
(j25l) as function of number of iterations for the Petviashvili method without accel¬ 
eration and for six acceleration techniqj^s with the best parameters mw and nw 
(according to Tables [5] and [H) . (b) Residual error (1251) as function of CPU time 
(in seconds) for six acceleration techniques with the best parameters mw and nw 
(^accordinff to Tables IHl and IHl) . ("cl Comnarison of residual error laEfainst CPU time) 

























nw 

AA-I(ntc) 

AA-II(nri;) 

1 

28(7.6888F - 14) 

30(3.6918A - 14) 

2 

21(3.2119A - 14) 

21(3.4864A - 14) 

3 

19(4.9065A - 14) 

20(4.6961A - 14) 

4 

19(2.5310F - 14) 

18(1.9550A - 14) 

5 

18(5.9489A - 14) 

18(4.2896A - 14) 

6 

18(2.0285F - 14) 

18(1.1674A - 14) 

7 

17(5.1055F - 14) 

17(5.5586A - 14) 

8 

17(2.8378A - 14) 

17(1.9837A - 14) 

9 

17(7.2231A - 14) 

16(7.1040A - 14) 

10 

17(2.8512A - 14) 

16(5.7507A - 14) 


Table 6 

Generalized solitary wave generation of (|37p . Number of iterations required by AA- 
I and AA-II as function of nw to achieve a residual error below TOL = 10“^^. The 
residual error at the last computed iterate is in parenthesis. Without acceleration, 
the Petviashvili method ([2]) , Q with 7 = 2 requires n = 53 iterations with a residual 
error 7.8594711 — 14. 

3.1.3 Numerical generation of periodic traveling waves of ^3^ 

The numerical generation of periodic traveling wave solutions of the BBM- 
BBM system fl38|) will complete the study about traveling wave generation 
of Boussinesq systems (127|) . Here the initial data are similar to those of the 
previous cases, although now Z = 16 is taken. Once system fIMll is solved, 
the application of the Petviashvili type method to fj35l) generates, for Ki = 
0.75,7^2 = 1 (taken as an example) the computed prohles shown in Figure 
[71(a)-(b). The periodic behaviour is also observed in the corresponding phase 
plots, shown in Figure [7|(c), (d), while the performance is illustrated in Figure 
[HI which corresponds to the behaviour of the residual error as function of the 
number of iterations. The method attains a residual error of 9.33536677 — 12 
in n = 572 iterations, showing the need of some acceleration technique. 

This slow performance is justified by the corresponding table of eigenvalues 
of the linearization operators, Table [7| in this case. We observe that besides 
eigenvalue one (associated to the translational invariance) the next largest 
in magnitude eigenvalue is close to one. (As in the generalized solitary wave 
generation the presence of conjugate complex eigenvalues, in this case with 
algebraic multiplicity above one, in the spectrum of the linearization matrices 
is noticed.) 

We now evaluate the application of VEM taking this example as illustration. 
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(a) (b) 




Fig. 7. Numerical generation of periodic traveling waves of (j38p . Approximate pro¬ 
files for Ki = 0.75, K 2 = 1. (a) rj profile; (b) u profile; (c) rj phase portrait, (d) u 
phase portrait. 



Fig. 8. Numerical generation of periodic traveling waves of (I38p . Approximate pro¬ 
files for Ki = 0.75, Ar 2 = 1. Residual error (j25l) vs number of iterations. 
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Iteration matrix S{gf,Uf) 

Iteration matrix F'{gf,Uf) 

2 . 000000 A -h 00 

l.OOOOOOA -b 00 

l.OOOOOOA -hOO 

-9.545242A - 01 

-9.545242A - 01 

-5.353103A - 01 - 6.459204A - Oli 

-5.353103A - 01 - 6.459204A - Olz 

-5.353103A - 01 - 6.459204A - Olz 

-5.353103A - 01 - 6.459204A - Oli 

-5.353103A - 01 -b 6.459204A - Oli 

-5.353103A - 01 -b 6.459204A - Olz 

-5.353103A - 01 -b 6.459204A - Oli 


Table 7 


Periodic traveling wave generation of (1381) with Ki = 0.75, 7^2 = 1- Six largest 
magnitude eigenvalues of the approximated iteration matrix S = L~^N'{r]f,Uf) 
(first column) and of the iteration matrix F'{r]f,Uf), generated by the Petviashvili 
method ([2]), (j4]) with 7 = 2, both evaluated at the last computed iterate {r]f,Uf). 

The standard comparison in performance is given in Table [3 In this case the 
tolerance for the residual error was set as TOL = 10“^^. Some conclusions 
from it are the following: 

(1) Better performance of polynomial methods compared to e-algorithms. 

(2) MPE and RRE are virtually equivalent, especially when k grows. There 
are more differences between VEA and TEA, but they decrease when k 
grows. 

(3) For polynomial methods, the best results are obtained for large k (around 
K = 9) while for e algorithms, it is better to take small k, (around mw = 
4, 5). This implies a similar length of each cycle (width of extrapolation). 

We now analyze the results corresponding to AAM by using Table |9l which 
evaluates the performance of AA-I and AA-II for the same example. Some 
conclusions from Table [9l 

(1) As in some previous cases the AAM (particularly AA-II) behave better 
than any VEM when measuring the performance in terms of the number 
of iterations. However, the polynomial methods MPE and RRE are more 
efficient in terms of the computational time, see Figure 13 

(2) The best results of AA-II are obtained with large values of nw. The 
method does not appear to be affected by ill-conditioning, contrary to 
AA-I, which becomes useless from nw = 7. 

3.2 Example 2. Localized ground state generation 


A second group of experiments illustrates the generation of localized ground 
states in nonlinear Schrodinger (NLS) type models with potentials. In partic- 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

1 

118 


278 

88 


(8.9607E - 12) 


(7.1073E - 12) 

(4.9649E - 12) 

2 

70 

81 

81 

81 


(9.2990E - 12) 

(6.8950E - 12) 

(6.8161E - 12) 

(6.8798E; - 12) 

3 

54 

53 

78 

64 


(9.4403E - 12) 

(7.2397E - 12) 

(4.2025E - 12) 

(9.4892E - 12) 

4 

47 

52 

55 

55 


(7.0369E - 12) 

(5.3720E - 12) 

(4.0063E - 12) 

(7.4587E - 12) 

5 

55 

46 

48 

103 


(4.0773E - 12) 

(8.2635E - 12) 

(9.7571E - 12) 

(6.2560E; - 12) 

6 

46 

44 

53 

79 


(9.9878E - 12) 

(5.3034E - 12) 

(5.7135E - 12) 

(5.4689E - 12) 

7 

45 

49 

51 

73 


(7.4645E - 12) 

(3.6062E - 12) 

(9.3703E - 12) 

(8.3118E - 12) 

8 

45 

46 

53 

69 


(9.1401E - 12) 

(4.3655E - 12) 

(9.0789E - 12) 

(4.4378E - 12) 

9 

41 

41 

58 

77 


(9.9555E - 12) 

(4.7100E - 12) 

(9.8141E - 12) 

(8.4583E - 12) 

10 

45 

45 

64 

64 


(3.9218E - 12) 

(4.4096E - 12) 

(6.2339E - 12) 

(5.0514E - 12) 


Table 8 

Periodic traveling wave generation of (I38p . Number of iterations required by MPE, 
RRE, VEA and TEA as function of k to achieve a residual error (1251) below TOL = 
10“^^. The residual error at the last computed iterate is in parenthesis. Without 
acceleration, the Petviashvili method ([2]), (|3|) with 7 = 2 requires n = 572 iterations 
with a residual error 9.3354E — 12. 

ular, the equation 


iut + dxxU + V {x)u + \u\^u = 0, (39) 

with potential V{x) = 6sech^(a:), is considered as an example, [i6ll80] . A lo¬ 
calized ground state solution of fl3^ has the form u{x,t) = e^^^U{x), where 
/i G M and U{x) is assumed to be real and localized {U —>■ 0, |a;| 00 ). 

Substitution into (l3^ leads to 
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nw 

AA-I(ntc) 

AA-II(nri;) 

1 

57(2.2757F - 12) 

66(8.3451U - 12) 

2 

84(4.7115F; - 12) 

40(6.4724U - 12) 

3 

81(6.2384F; - 12) 

39(3.1769U - 12) 

4 

81(1.4968F - 12) 

36(3.2197U - 12) 

5 

48(9.8634F; - 12) 

36(4.5059U - 12) 

6 

48(3.2036F - 12) 

49(3.7796U - 12) 

7 

Ill-conditioned 

38(2.4939U - 12) 

8 

Ill-conditioned 

36(7.2839U - 12) 

9 

Ill-conditioned 

35(8.3908U - 12) 

10 

Ill-conditioned 

37(2.6179U - 12) 


Table 9 

Periodic traveling wave generation of (I38p . Number of iterations required by AA-I 
and AA-II as function of nw to achieve a residual error (1251) below TOL = 10“^^. 
The residual error at the last computed iterate is in parenthesis. Without accelera¬ 
tion, the Petviashvili method ([2|), (jH) with 7 = 2 requires n = 572 iterations with a 
residual error 9.3354711 — 12. 



Fig. 9. Numerical generation of periodic traveling waves of pSI) . Residual error (|25p 
as function of: (a) number of iterations and (b) CPU time (in seconds) for MPE(9) 
and AA-II(9). 

U''{x) + V{x)U{x) - iiU{x) + \U{x)\^U{x) = 0. (40) 

A discretization of (1401) based on a Fourier collocation method on a sufficiently 
long interval (—/, 1) requires in this case the resolution of a system of the form 
© for the approximations Uh of U at the grid points Xj = —I + jh, h = 
2l/m,j = 0,..., m — 1, with 
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L = D‘^ + diagiV) - N{Uh) = -Uh.\ 


where D is the pseudospectral differentiation matrix, diag(id) is the diagonal 
matrix with elements Vj = V(xj),j = 0,... ,m — 1 and Im is the m x m 
identity matrix. The nonlinearity N is homogeneons with degree three, where, 
as usnal, the dot stands for the Hadamard product. The discussion below is 




(a) (b) 



(c) (d) 


Fig. 10. Numerical generation of localized ground states of (j40p . Approximate asym¬ 
metric profile, (a) p = 1.3; (b) p = 3.3; (c) p = 6.3; (d) p = 8.3. The amplitude of 
the waves increases with p, while the shape is narrower. 

focused on the ground state numerical generation for several values of p, which 
provide different challenges to the iteration. For each considered value of p, 
the performance of both families of acceleration techniques has been checked. 

The first results concern the numerical generation of an asymmetric solution 
of flTUD for p = 1.3 iFigure [TUfa)!. Figure [TTl compares the performance of 
the acceleration techniques in terms of the number of iterations required to 
reduce the residual error fl2^ below TOL = 10“^^ and as function of the ex¬ 
trapolation width parameters k and nw. In the case of VEM, Figure m^), 
all the techniques considered are comparable and the differences are not very 
large; MPE with k = 7 gives the minimum number of iterations. (The values 
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K = 8, 9 also lead to the same number of iterations, but the computational 
effort in CPU time is higher.) The AAM, Figure fTTT bh are competitive with 
VEM for small values of nw {nw = 1,2). As nw grows, the number of iter¬ 
ations increases (in opposite way to the behaviour of VEM with respect to 
k) and the computation of the coefficients in the minimization problem be¬ 
comes ill-conditioned. The comparison between the most efficient method of 
each family (MPE(7) and AA-II(2) respectively) is displayed in Figure [121 It 



Fig. 11. Numerical generation of asymmetric ground state of (I40p with p = 1.3. 
Number of iterations required to reduce the residual error (I25p below TOL = 10“^^ 
and as function of the extrapolation width parameters k and nw. (a) VEM; (b) 
AAM. 

shows the residual error as function of the number of iterations (a) and the 
CPU time (b) for the Petviashvili method without acceleration (solid line) 
and accelerated with MPE(7) and AA-II(2). In both hgures, the improvement 
in the performance with respect to the Petviashvili method provided by the 
two acceleration techniques is observed, with the best results corresponding to 
MPE (and, in general VEM against AAM). Table [TOl conhrms the convergence 
of the Petviashvili method. It displays the six largest magnitude eigenvalues 
of the corresponding iteration matrix of the classical hxed-point algorithm 
S = L~^N'{Uf), and of the Petviashvili method ([3]), (01) for two values of p. 
Since an analytical expression for the exact prohle is not known, the matrices 
have been evaluated at the last computed iterate given by MPE(7). In the 
case of S (hrst column), the dominant eigenvalue corresponds to the degree 
of homogeneity p = 3, with the rest of the eigenvalues below one. The hlter 
action of the stabilizing factor is observed in the second column. The degree 
p = 3 has been subtituted by zero (the optimal <? = 7(1 — p) = —p has been 
taken) and the rest of the spectrum is preserved. This implies that for p = 1.3, 
the spectral radius of F'{u*) is below one (second column) and this leads to 
the (local) convergence of the method. 

For other values of p, some differences are observed. When p = 3.3 the numer- 
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(a) 



Fig. 12. Numerical generation of assymetric ground state of (I40p with fj, = 1.3. Resid¬ 
ual error (I25p as function of the number of iterations (a) and CPU time in seconds 
(b) for the Petviashvili method without acceleration (solid line) and accelerated 
with MPE(7) (circle symbols) and AA-II(2) (plus symbols). 

ical generation of an asymmetric solntion of (HOI) (see Figure fTOf b)) with the 
Petviashvili method withont acceleration is not possible in general. Table [10] 
(third column) shows the presence of an additional eigenvalue with magnitude 
above one in the iteration matrix S of the classical hxed point algorithm. As 
part of the spectrum different from the degree of homogeneity p = 3, this 
eigenvalue also appears in the spectrum of the iteration matrix of the Petvi¬ 
ashvili method (fourth column) making thus the convergence fail. Here the use 
of the acceleration techniques corrects this behaviour, leading to convergence. 
(Both iteration matrices are in fact evaluated at the approximate prohle dis¬ 
played in Figure fTOf b).) In this case (see Figure [T3f a)) MPE and RRE have 


37 
















p = 1.3 

p = 1.3 

II 

CO 

CO 

p = 3.3 

eigs(S') 

eigs{F' (u*)) 

eigs(S’) 

eigs{F' (u*)) 

2.999999E-F00 

2.886842E-01 

-6.328271E-F00 

-6.328271E-F00 

2.886842E-01 

-1.858331E-01 

3.000000E-h00 

8.594730E-01 

-1.858331E-01 

1.419117E-01 

8.594730E-01 

5.552068E-01 

1.419117E-01 

7.522396E-02 

5.552068E-01 

2.978699E-01 

7.522396E-02 

5.527593E-02 

2.978699E-01 

2.360730E-01 

5.527593E-02 

3.629863E-02 

2.360730E-01 

1.552434E-01 


Table 10 

Numerical generation of asymmetric profile of (|40p with /i = 1.3 and ^ = 3.3. Six 
largest magnitude eigenvalues of the approximated iteration matrix of the classical 
fixed-point method S = L~^N' {Uf) (left) and of the Petviashvili method, evaluated 
at the last computed iterate Uj obtained with MPE(7). 

virtually the same performance while the e-algorithms start reducing their 
efficiency. (In this case, TEA does not always work in a reliable way and is 
not competitive against the other VEM.) As far as the AAM are concerned. 
Figure [T3f b). both improve the performance in a similar, relevant way. They 
are comparable with VEM in number of iterations (Figure fTTf a)) and behave 
better when measuring the computational time (Figure [TU(b)). 




(a) (b) 

Fig. 13. Numerical generation of asymmetric ground state of (l4Up with p = 3.3. 
Number of iterations required to reduce the residual error (j25p below TOL = 10“^^ 
and as function of the extrapolation width parameters n and nw. (a) VEM; (b) 
AAM. 


It may be worth considering the case p = 6.3 because of some relevant points. 
The hrst one is the generation of the asymmetric prohle. Figure flOf c). which 
in general is not possible with the Petviashvili method without acceleration. 
The situation is similar to that of the previous case p = 3.3 and it is shown 
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Fig. 14. Numerical generation of asymmetric ground state of (14011 with /x = 3.3. 
Residual error (j25p as function of the number of iterations (a) and CPU time in 
seconds (b) for the Petviashvili method without acceleration (solid line) and accel¬ 
erated with MPE(8) (circle symbols) and AA-II(5) (plus symbols). 

in Table [H] (first and second columns). In this case, the best results of the 
acceleration are given by MPE and AA-I (Figure [T5ll . The loss of performance 
of the e-algorithms and the improvement of AAM, observed in the previous 
experiments, are conhrmed here and in the experiments for /i = 8.3 (Figures 
EH and fT^ . The comparison between MPE and AA-I, see Figures fTHf a). (b), 
reveals, ikn the authors’ opinion, a similar performance. 



(a) (b) 

Fig. 15. Numerical generation of asymmetric ground state of (14011 with p = 6.3. 
Number of iterations required to reduce the residual error (I25p below TOL = 10“^^ 
and as function of the extrapolation width parameters k and nw. (a) VEM; (b) 
AAM. 


The second question with regard to the case /x = 6.3 concerns the behaviour 
of the Petviashvili method without acceleration. In this case, the method is 
convergent, but to a symmetric localized wave, see Figure [191 This can be 
explained by the first two columns of Table [TT] and by Table EH Note that, as 
mentioned before, for the asymmetric solution, the Petviashvili method can- 
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Fig. 16. Numerical generation of asymmetric ground state of (I40p with /x = 6.3. 
Residual error (j25p as function of the number of iterations (a) and CPU time in 
seconds (b) for the Petviashvili method without acceleration (solid line) and accel¬ 
erated with MPE(IO) (circle symbols) and AA-I(5) (triangle symbols). 


/X = 6.3 

/X = 6.3 

II 

OO 

CO 

CO 

00 

II 

eigs(S’) 

eigs(F'(M*)) 

eigs(S’) 

eigs(F'(xi*)) 

5.095370E-h00 

5.096207E-F00 

3.962824E-h00 

3.962824E-F00 

3.000000E-h00 

9.672929E-01 

2.999999E-h00 

9.807797E-01 

9.672929E-01 

7.506018E-01 

9.807797E-01 

8.081404E-01 

7.506018E-01 

4.078905E-01 

8.081404E-01 

4.459845E-01 

4.078905E-01 

3.472429E-01 

4.459845E-01 

4.030040E-01 

3.472429E-01 

2.032986E-01 

4.030040E-01 

1.929797E-01 


Table 11 

Nnmerical generation of asymmetric profile of (I40p with ^ = 6.3 and /U = 8.3. Six 
largest magnitnde eigenvalues of the approximated iteration matrix of the classical 
fixed-point method S = L~^N'{Uf) (left) and of the Petviashvili method (l2|), ([3]), 
evaluated at the last computed iterate Uf obtained with MPE(IO). 

not be convergent. However, according to the information provided by Table 
m this is locally convergent to the symmetric solution. (In this case, the spec¬ 
tral radius of F'{u*) is below one.) This profile can be indeed approximated 
by using acceleration techniques (and with the corresponding computational 
saving) but starting from a different initial iteration. 

Finally, the case /i = 8.3 is also analyzed, see Figure ITOl d). The main reason 
we find to emphasize this case is to confirm the conclusions obtained from the 
experiments with the previous values of /i: 

• Among the VEM, the polynomial methods give a better performance, while 
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eigs(S') 

e\gs{F'{u*)) 

3.000000E-F00 

6.098684E-01 

6.098684E-01 

2.696853E-01 

2.696853E-01 

1.518421E-01 

1.518421E-01 

1.039553E-01 

1.039553E-01 

6.046185E-02 

6.046185E-02 

5.492737E-02 


Table 12 

Numerical generation of symmetric profile of (I40p with ^ = 6.3. Six largest mag¬ 
nitude eigenvalues of the approximated iteration matrix of the classical fixed-point 
method S = L~^N'{Uf) (left) and of the Petviashvili method dS]), ([3]), evaluated at 
the last computed iterate Uf obtained with Petviashvili method (1^ . ll4l) . 


the e-algorithms become less efficient as /i increases. As observed in Figure 
[Tnl the larger /i the larger and narrower the asymmetric prohle is. The com¬ 
putation becomes harder as is noticed by comparing the iterations required 
by the methods in Figures [11], [131 [El cind [TT] One can also note the in¬ 
crement of the magnitude of the eigenvalues of the corresponding iteration 
matrices of the Petviashvili method in Tables [10] and [HI Therefore, under 
more demanding conditions, the polynomial methods give a better answer 
than the e-algorithms. 

• Contrary to the e-algorithms, whose performance gets worse as /i increases, 
the AAM improve their behaviour up to being comparable with polynomial 
methods (cf. the periodic traveling wave generation in Section [3]). Further¬ 
more, this is obtained with small values of the parameter nw, thus avoiding 
ill-conditioned problems. 


3.3 Example 3. Solitary wave solutions of the Benjamin equation 


An additional application of acceleration techniques concerns the oscillatory 
character of the wave to be numerically generated. This property has shown 
to make influence on the performance of the iteration with eventual loss of 
convergence in some cases, [3IE2]. Presented here is the use of acceleration as 
an alternative to overcome this difficulty and this will be illustrated with the 
numerical generation of solitary waves in one- and two-dimensional versions 
of the Benjamin equation. 
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(a) (b) 

Fig. 17. Numerical generation of asymmetric ground state of (I40p with /x = 8.3. 
Number of iterations reqnired to reduce the residual error (j25h below TOL = 10“^^ 
and as function of the extrapolation width parameters k and nw. (a) VEM; (b) 
AAM. 



(a) (b) 


Fig. 18. Numerical generation of asymmetric gronnd state of (I40p with /x = 8.3. 
Residual error (j25p as function of the number of iterations (a) and CPU time in 
seconds (b) for the Petviashvili method without acceleration (solid line) and accel¬ 
erated with MPE(7) (circle symbols) and AA-II(3) (plus symbols). 

3.3.1 One-dimensional Benjamin equation 

A first example of the situation described above is given by the solitary wave 
solutions of the Benjamin equation, [7j 

tit -\- OCVjx [dtltlx 'yB.Vjxx dUxxx O7 

where u = u{x, t),x E > 0, a, (3,^,6 are positive constants, and Ti denotes 
the Hilbert transform dehned on the real line as 

nf{x) :=-p.v. f (42) 

TT J -00 X — y 
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Fig. 19. Numerical generation of symmetric ground state of (j40h with ^ = 6.3. 
Approximate profile with Petviashvili method without acceleration. 

or through its Fourier transform as 


Hfik) = —isign{k)f{k), k eR. 


Equation fHT]) is a model for the propagation of internal waves along the 
interface of a two-layer fluid system and where gravity and surface tension 
effects are not negligible. It includes the limiting cases of negligible surface 
tension (5 = 0 or Benjamin-Ono equation) and a limit of a model with very 
thin upper fluid (7 = 0 or KdV equation). Solitary-wave solutions of fHTll 
with speed c* > 0 are determined by prohles u{x, t) = ip{x — Cgt), Cg > 0, such 
that ip and its derivatives tend to zero as X = x — Cgt approaches ±cxo and 
satisfying 



(43) 


where ' = d/dX. Albert et ah, [2] established a complete theory of existence 
and orbital stability of solitary waves of fHT]) for small 7, while Benjamin, [8], 
derived the oscillating behaviour of the waves, with the number of oscillations 
increasing as 7 approaches 7* = 5{a — c^), along with the asymptotic decay, 

as |X| —)■ 00, like 1/X^. 

Except in the limiting cases, solitary wave solutions are not analytically known. 
A standard way to generate solitary wave prohles numerically consists of con¬ 
sidering in the Fourier space 



(44) 
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(where ^{k) is the Fourier transform of ip), discretizing with periodic 
boundary conditions on a sufficiently long interval (—/, /) and the use of dis¬ 
crete Fourier transform (DFT) 


{-Cs +a-'y\k\+5k‘^)ip^i. + ^{^^ *ip^^ =0, (45) 

for A; = — Y,..., Y — 1, where tp^ is a trigonometric polynomial of degree N 
which approximates ip and denotes its Fourier coefficient. Then (1451) is 
numerically solved by incremental continuation with respect to 7 from 7 = 0 
(which corresponds to KdV equation and for which solitary wave prohles are 
analytically known) and a nonlinear iteratively solver for each value of the 
homotopic path with respect to 7 . For a more detailed description of the 
incremental continuation method and the performance of several nonlinear 
iterative solvers see mn- The experiments performed there reveal that the 
oscillatory behaviour of the wave increases the difficulty of its computation, 
even using numerical continuation. Our aim here is giving a computational 
alternative, based on the use of acceleration techniques. 

To this end, we £x a speed Cg = 0.75, the parameters a = (3 = 6 = 1 and 
generate numerically a solitary wave solution of fl43|) by combining the Petvi- 
ashvili method, standing for the family of iterative method ([ 2 ]), along with the 
acceleration techniques considered in previous examples. We will take four 
values of 7 , namely 0.9, 0.99, 0.999, 0.9999 (which, for the considered values of 
the parameters, are close to 7 *, equals 1 in our example), correspond to a more 
and more oscillating prohle (with smaller and smaller amplitude, see Figures 
[20]f a)-fd): the computational window is [—512, 512] with N = 4096 collocation 
points) and for which the Petviashvili method with numerical continuation re¬ 
quires a long computation to converge or directly does not work. In all the 
experiments the initial iteration is the (analytically known) solitary wave pro¬ 
hle corresponding to 7 = 0 (KdV equation). As in the previous examples, we 
hrst estimate the performance of the acceleration techniques by comparing 
the number of iterations required by each of them to achieve a residual error 
(12 5 p less than a tolerance TOL = 10“^^. For the case of the VEM and the 
four values of 7 considered, this information is given in Tables [T41fTBl All the 
methods achieve convergence in the four cases. (The Petviashvili method with 
continuation is not able to converge for the last two values of 7 and for the hrst 
two values the number of iterations required is prohibitive: for example, just 
going from 7 = 0.98 to 7 = 0.99 the method requires 266 iterations to have a 
residual error of size 9.063477 — 14; the continuation process from the initial 
7 = 0 , where our computations start, requires a total number of iterations of 
about 4470.) As expected the ehort of VEM in number of iterations increases 
with 7 , that is, with the oscillating character of the prohle, see Figure l?IT a). 
Among them and except some particular cases (for example, when TEA is ap- 
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0.06 



Fig. 20. Approximate solitary-wave solutions of (|43p given by MPE method with 
K, = 7 and Cg = 0.75, a = /3 = 5 = 1. (a)-(d) correspond, respectively, to 
7 = 0.9,0.99,0.999,0.9999. 


plied with 7 = 0.9999 and k = 6) the polynomial methods are more efficient 
than e-algorithms when 7 increases, althongh the difference is shorter than 
that was obtained in the examples of Section 13.11 It is remarkable that in the 
case of a solitary wave prohle with a small nnmber of oscillations (for example, 
when 7 = 0.9) the performance of the methods is virtnally the same: after one 
or two cycles, the improvement of the acceleration techniqne is good enongh to 
not needing to complete the next cycle in order to achieve the tolerance for the 
residnal error. This is particularly emphasized in the case of the e-algorithms, 
where the cycle is longer {2k against k + 1 for the case of the polynomial 
methods). The results from Tables fT3llT6] are in contrast with those from Ta¬ 
bles [T71I201 that correspond to the AAM. The main conclusion here is that 
these methods are strongly affected by the oscillating character of the prohles 
and, compared to VEM, do not seem to be recommendable for this sort of 
computations, at least without a suitable choice of preconditioning. (Some of 
it was suggested by the previous experiments concerning the generalized soli¬ 
tary waves of some Boussinesq systems, see Section O) Tables 1171120] show that 
as 7 increases, ill-conditioning of the corresponding least-squares problem is 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

2 

35 

35 

33 

36 


(7.0418E - 14) 

(5.3048E - 14) 

(7.3270E - 14) 

(7.5208E - 14) 

3 

26 

28 

28 

25 


(8.5843E - 14) 

(3.7474E - 14) 

(4.3193E - 14) 

(3.5514E - 14) 

4 

25 

27 

23 

21 


(8.7162E - 14) 

(4.0422E - 14) 

(2.6305E - 14) 

(2.7078E - 14) 

5 

22 

22 

25 

25 


(2.3778E - 14) 

(2.4540E - 14) 

(2.4551E - 14) 

(2.6612E - 14) 

6 

19 

19 

24 

24 


(6.8569E - 14) 

(6.3344E - 14) 

(8.5839E - 14) 

(8.8986E - 14) 

7 

19 

19 

21 

24 


(4.8718E - 15) 

(4.4709E - 15) 

(5.1694E - 14) 

(7.2978E - 14) 

8 

21 

21 

19 

19 


(2.4289E - 15) 

(2.3872E - 15) 

(1.8080E - 14) 

(6.5383E - 14) 

9 

23 

23 

21 

21 


(2.2138E - 15) 

(2.5529E - 15) 

(3.5794E - 15) 

(3.3103E - 15) 

10 

25 

25 

23 

23 


(2.3084E - 15) 

(2.2764E - 15) 

(4.3858E - 15) 

(3.3545E - 15) 


Table 13 

Solitary wave generation of (|43p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error (I25|) below TOL = 
10“^^. The residual error at the last computed iterate is in parenthesis; Cg = 0.75, 
a = f3 = S=l, j = 0.9. 

observed from even moderate values of nw in the case of AA-I (affecting the 
stability of the method, which is not able to converge or requires a great effort 
in number of iterations) while AA-II is not so affected. However, when AAM 
work, they exhibit a competitive performance, as shown in Figure [2TT bl when 
comparing with Figure |2T](a). (Our implementation follows that of described 
in na, which uses the unconstrained form of the least-squares problem and 
that also was suggested by some other authors, [M]. For the numerical res¬ 
olution we have used some other alternatives, with QR decomposition with 
pivoting, [75], and SVD, [Mj. The results of Tables ITTII^ correspond to the 
hrst implementation, while the second one overcomes ill-conditioning in some 
more cases of AA-I, but at the cost of an important increase of the iterations.) 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

2 

77 

67 

147 

93 


(4.3301E - 14) 

(5.8218E - 14) 

(1.8199E - 14) 

(8.9274E - 14) 

3 

56 

53 

47 

51 


(2.9502E - 14) 

(4.4310E - 14) 

(2.1373E - 14) 

(6.5537E - 14) 

4 

43 

49 

48 

41 


(4.7982E - 14) 

(1.5692E - 14) 

(7.4855E - 14) 

(2.4388E - 14) 

5 

38 

36 

39 

37 


(3.8410E - 14) 

(2.5110E - 14) 

(4.0062E - 14) 

(3.4320E - 14) 

6 

33 

33 

43 

43 


(3.6079E - 14) 

(2.3031E - 14) 

(1.6073E - 14) 

(1.6545E - 14) 

7 

30 

30 

35 

49 


(6.3053E - 14) 

(3.2262E - 14) 

(3.3824E - 14) 

(1.8129E - 15) 

8 

38 

31 

40 

37 


(9.5443E - 14) 

(6.3679E - 14) 

(3.2121E - 14) 

(1.4194E - 14) 

9 

34 

34 

41 

41 


(1.8749E - 15) 

(2.1169E - 15) 

(4.1825E - 15) 

(1.8746E - 15) 

10 

37 

37 

45 

45 


(2.2503E - 15) 

(1.7597E - 15) 

(4.7245E - 15) 

(1.8872E - 15) 


Table 14 

Solitary wave generation of (|43p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error (I25|) below TOL = 
10“^^. The residual error at the last computed iterate is in parenthesis; Cg = 0.75, 
a = f3 = S=l, j = 0.99. 

3.3.2 Lump solitary waves of 2D Benjamin equation 

In order to finish off this example we stndy the performance of the acceleration 
techniqnes when generating nnmerically Inmp solitary wave solntions of the 
2D Benjamin equation ^311441145] 

{rjt + = 0, (46) 

where a, /3,6 > 0 and "H is the Hilbert transform (I42|) with respect to x. In 
fl46|) . as in the one-dimensional case, rj = ri{x,z,t) stands for the interfacial 
deviation wave between two ideal fluids with a bounded upper layer and the 
heavier one with inhnite depth, and under the presence of interfacial tension. 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

3 

63 

60 

56 

65 


(2.2431E - 14) 

(3.9572E - 14) 

(9.8315E - 14) 

(6.6026E - 14) 

4 

79 

49 

77 

70 


(4.9749E - 14) 

(4.1392E - 14) 

(6.6180E - 14) 

(9.3318E - 14) 

5 

43 

43 

54 

101 


(5.8270E - 14) 

(1.0725E - 14) 

(5.9093E - 14) 

(4.4323E; - 14) 

6 

45 

45 

71 

71 


(9.4355E - 14) 

(3.5010E - 14) 

(6.8677E - 14) 

(1.9108E - 14) 

7 

37 

37 

65 

65 


(8.3020E - 14) 

(6.3139E - 14) 

(1.1800E - 14) 

(2.2238E - 15) 

8 

51 

41 

73 

73 


(4.8702E - 15) 

(1.3712E - 14) 

(8.5215E - 15) 

(3.7373E - 15) 

9 

45 

45 

81 

81 


(2.8855E - 15) 

(2.7603E - 15) 

(3.4268E - 15) 

(2.6722E - 15) 

10 

51 

51 

111 

69 


(3.3110E - 14) 

(3.2354E - 14) 

(2.0694E - 15) 

(4.0083E - 14) 


Table 15 

Solitary wave generation of (I43p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error (j25h below TOL = 
10“^^. The residual error at the last computed iterate is in parenthesis; Cg = 0.75, 
a = j3 = 6=l,j = 0.999. 




(a) (b) 

Eig. 21. Solitary wave generation of (|43]) . Cg = 0.75, a = ^ = 6 = 1. Number of itera¬ 
tions required to achieve a residual error (in Euclidean norm) below TOL = IE —13 
as function of 7. (a) MPE, RRE, VEA and TEA. (b) AA-I and AA-II. 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

3 

60 

56 

65 

49 


(4.0201E - 14) 

(2.3990E - 14) 

(1.2911E - 14) 

(1.4934E - 14) 

4 

55 

61 

203 

51 


(2.4890E - 14) 

(4.0300E - 14) 

(6.8061E - 14) 

(1.4970E - 14) 

5 

45 

42 

52 

73 


(1.4175E - 14) 

(6.6916E - 14) 

(5.7187E - 14) 

(1.2493E; - 15) 

6 

43 

41 

101 

43 


(3.6348E - 14) 

(1.5143E - 14) 

(7.7548E - 14) 

(4.9096E - 14) 

7 

46 

46 

65 

65 


(4.7457E - 15) 

(2.3306E - 15) 

(9.9131E - 14) 

(8.3133E - 16) 

8 

51 

51 

109 

73 


(1.9157E - 15) 

(7.6397E - 16) 

(3.0772E - 15) 

(1.3583E; - 15) 

9 

56 

56 

101 

81 


(6.7580E - 16) 

(7.1896E - 16) 

(7.8767E - 16) 

(2.4959E - 15) 

10 

61 

61 

199 

89 


(7.8194E - 16) 

(9.0591E - 16) 

(7.7246E - 16) 

(6.4547E - 16) 


Table 16 

Solitary wave generation of (I43p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error (l25p below TOL = 10“^^. 
The residual error at the last computed iterate is in parenthesis; Cg = 0.75, a = (3 = 
5 = 1, 7 = 0.9999. 


nw 

AA-I(nrc) 

AA-II(nrc) 

1 

20(5.5496E - 14) 

25(4.1169E - 14) 

2 

17(1.9530E - 15) 

21(1.2546E - 14) 

3 

15(4.4608E - 15) 

20(4.8161E - 14) 

4 

15(1.3401E - 14) 

20(8.3292E - 14) 

5 

14(5.1185E - 14) 

21(1.4961E - 14) 

6 

14(1.5211E - 14) 

21(1.4958E - 14) 


Table 17 

Solitary wave generation of ()43p . Number of iterations required by AA-I and AA-II 
as function of nw to achieve a residual error (I25p below TOL = 10“^^. The residual 
error at the last computed iterate is in parenthesis. Cg = 0.75, a = /3 = b = 1, 
7 = 0.9. 
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nw 

AA-I(ntc) 

AA-II(nrc) 

1 

58(8.4209F - 14) 

38(6.6459.F - 14) 

2 

44(3.1097F; - 14) 

44(5.6598.F - 14) 

3 

38(1.0963F; - 14) 

36(9.8103.F - 14) 

4 

53(2.9810F - 14) 

47(3.1570.F - 14) 

5 


57(4.4383.F - 14) 

6 


47(5.3103.F - 14) 

7 


53(7.2250.F - 14) 

8 


64(3.6054.F - 14) 


Table 18 

Solitary wave generation of (I43p . Number of iterations required by AA-I and AA-II 
as function of nw to achieve a residual error ()25p below TOL = 10“^^. The residual 
error at the last computed iterate is in parenthesis. Cg = 0.75, a = /3 = 6 = 1, 
7 = 0.99. 


nw 

AA-I(ntc) 

AA-II(ni/;) 

1 

63(8.2392F - 14) 

77(9.2229.F - 14) 

2 

57(9.5654F; - 14) 

83(8.6540.F - 14) 

3 

225(4.7416.F - 14) 

49(7.4595.F - 14) 

4 

33(2.0366F - 14) 

91(6.5926.F - 14) 

5 


64(1.1616.F - 14) 

6 


70(5.9306.F - 14) 


Table 19 

Solitary wave generation of (143 p . Number of iterations required by AA-I and AA-II 
as function of nw to achieve a residual error (I25p below TOL = 10“^^. The residual 
error at the last computed iterate is in parenthesis. Cg = 0.75, a = (3 = <5 = 1, 
7 = 0.999. 

The two-dimensional version incorporates weak transverse variations. For the 
experiments below we will consider a normalized version of (06]), |13] 


(hi + ijf)x - 2VU{r]^^) + rjxxx)^ - Vzz = 0, (47) 

where F > 0. (The case F = 0 corresponds to the KP-I equation, [12].) For 
localized solutions, the zero total mass 


ri{x, z, t)dx = 0, 


(48) 
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nw 

AA-I(nt(;) 

AA-II (ntc) 

1 

800(2.7029L - 14) 

84(9.3995L - 15) 

2 

61(3.1858L - 14) 

193(3.5965L - 14) 

3 

99(2.2611L - 14) 

98(1.3640L - 14) 

4 

144(8.0198L - 14) 

82(4.5365L - 14) 

5 


94(1.8114L - 14) 

6 


66(8.1966L - 14) 


Table 20 

Solitary wave generation of (143 p . Nnmber of iterations required by AA-I and AA-II 
as function of nw to achieve a residual error (I25p below TOL = 10“^^. The residual 
error at the last computed iterate is in parenthesis. Cg = 0.75, a = P = S = 1, 
7 = 0.9999. 

is also assumed. Lump solitary wave solutions of (HTl) . (HHl) are solutions of the 
form 7]{x,z,t) = ri{X,Z),X = x — Cst,Z = 2 ; for some Cg > 0. Substitution 
into glD leads to 


(-Cgi] + rf - 2VH{r]x) + Vxx)- Vzz = 0, (49) 

As shown in [15], the value L = 1 marks a bifurcation point as for the type 
of lump solutions of (1471) between lumps of KP-I type and of wavepacket 
type. This implies in particular that as T < 1 approaches one the lump wave 
increases the oscillations. 

The numerical procedure used in [43B5] to generate approximate lump waves 
combines numerical continuation in T, pseudospectral approximation to (1491) 
(where constraint (j48l) is imposed) and Newton’s method for the resolution of 
the corresponding system of equations in each step of the T-homotopic path. 
The use of the Petviashvili methods (instead of Newton’s) was suggested in [3]. 
(For the use of the Petviashvili method in the generation of two-dimensional 
solitary waves see e. g. [nzu.) As in the one-dimensional case, the compu¬ 
tation of approximate lump prohles comes up two main difficulties: the use 
of numerical continuation and the oscillating behaviour of the lump. These 
problems can be overcome with the use of acceleration techniques, especially 
VEM. 

In order to illustrate this we will take Cg = 1 and generate approximate lump 
solitary waves for P = 0.99,0.999,0.9999. As described in [5], the periodic 
problem on a square [—La,,L 3 ;] x [—Lz,Lz] of (147)) is discretized by using a 
Fourier collocation method, generating approximations {r]h)ij to the lump 
prohle ri{xi, Zj) at the collocation points x* = —La, + ihx, Zj = —L^ + jh^, hg. = 
2Lx/Nx, hz = 2Lz/Nz, i = I,..., Nz;,j = 1,..., Nz. The system for the discrete 
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Fourier coefficients of the approximation is of the form 


{kl{cs + 2V%\ + kl))^h{kx, k,) = kl (rif^ {k^,, k,), (50) 

for kx = —Nx/2,..., Nx/2,kz = —Nz/2 ,..., iV^/2 and where fjh{kx, k^) stands 
for discrete /c^)-Fourier component of rjh- The zero total mass condition 
(HHi) is imposed as 


%(0,0) = 0. (51) 

When (1^ is included into fl50H . the resulting system for the rest of Fourier 
components is nonsingular and it is iteratively solved, for hxed F, by using: 


(i) The Petviashvili method with numerical continuation from the initial iter¬ 
ation given by the exact prohle for F = 0 


rio{x,z) 


I2cs 


3 + cl z'^ — CsX^ 
(3 -|- + clz‘^y 


(52) 


(ii) The Petviashvili method without numerical continuation but accelerated 
with the six techniques MPE, RRE, TEA, VEA, AA-I and AA-II and the 
same initial iteration ([52]). 


The experiments below follow a similar design to that of the one dimensional 
case. We have taken = 1024 with = 256 and a toler¬ 

ance of TOL = 10“^ for the control of the iteration. As before, the number 
of iterations shown in the numerical results correspond to the total account, 
including the iterations of each cycle. From this value, one can obtain the it¬ 
erations exclusively due to the corresponding acceleration. We think that this 
way of counting the iterations makes the comparison with the results without 
acceleration more realistic. We also remark that the use of the same initial 




iteration (|52|) is against the alternative technique with acceleration methods 
since, according to the form of the resulting waves, the initial prohle is not 
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Fig. 22. Solitary wave generation of (|47p with Petviashvili method, accelerated with 
MPE. The approximated profiles correspond to P = 0.99,0.999,0.9999 (left). On the 
right, the corresponding X and Z cross sections are shown (solid and dashed-dotted 
lines, respectively) . 

close. This should be observed in the behaviour of the residual error with re¬ 
spect to the number of iterations: the main effort is at the beginning; once 
the error is small enough, all the techniques accelerate the convergence in a 
more important way. Here three values of T = 0.99, 0.999, 0.9999 are consid¬ 
ered. The corresponding approximate waves (conhrming the highly oscillatory 
behaviour) are computed with the acceleration procedure of (ii) and can be 
observed in Figure [22l The hrst procedure in (i), based on continuation with 
respect to T, is totally inefficient for the the hrst value and does not work for 
the other two. The performance of the acceleration techniques is compared in 
Figure [23] and Tables l2Tli2ll 

The comparison of the techniques in this case conhrms the conclusions ob¬ 
tained in the one-dimensional version, namely: 

• The best performance is given by the polynomial methods (MPE in this 
case). 

• The e-algorithms, although less efficient, are also competitive (contrary to 
what was observed in some previous examples). 

• AAM only work correctly up to a moderate value of F < 1. When F ap¬ 
proaches one they cannot get the performance of VEM or directly fail. 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

2 

38 

41 

63 

49 


(7.5320E - 09) 

(7.7910E - 09) 

(5.8233E - 09) 

(8.1460E - 09) 

3 

32 

75 

49 

56 


(6.6515E - 09) 

(3.3504E - 09) 

(9.0167E - 09) 

(9.4419E - 09) 

4 

32 

37 

33 

38 


(5.0660E - 09) 

(3.8941E - 09) 

(9.4694E - 09) 

(8.9305E - 09) 

5 

26 

30 

35 

34 


(4.9957E - 09) 

(8.4036E - 09) 

(1.3179E-09) 

(1.1643E -09) 

6 

29 

30 

30 

27 


(4.5762E - 09) 

(1.4663E -09) 

(9.7151E - 09) 

(2.8493E - 09) 

7 

25 

33 

31 

31 


(2.8813E - 09) 

(5.5247E - 09) 

(1.8370E - 09) 

(6.0375E - 09) 

8 

21 

29 

35 

35 


(6.8132E - 09) 

(3.7396E - 09) 

(6.7544E - 09) 

(2.3207E - 09) 

9 

22 

25 

35 

40 


(7.0811E - 09) 

(6.1487E - 09) 

(7.6658E - 09) 

(5.3006E - 09) 

10 

23 

25 

39 

43 


(9.6176E - 10) 

(7.8664E - 09) 

(9.0031E - 09) 

(2.7713E - 09) 


Table 21 

Solitary wave generation of (|47p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error below TOL = 10“®. The 
residual error ([25|) at the last computed iterate is in parenthesis; = 1, T = 0.99. 

4 Acceleration techniques with extended Petviashvili type meth¬ 
ods 

One of the drawbacks of the Petviashvili type methods ([2]) in traveling wave 
generation is their limitation to some specihc problems, namely those with ho¬ 
mogeneous nonlinearities. When the nonlinear term is not homogeneous but 
a combination of homogeneous functions of different degree, these methods 
can be extended by adapting the stabilizing function s to each homogeneous 
part. This leads to the so-called e-Petviashvili type methods, derived in [1]. In 
this section and in order to improve the traveling wave generation for prob¬ 
lems with this type of nonlinearities, we will apply the acceleration techniques 
to the e-Petviashvili method as initial iterative procedure. This will be illus- 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

3 

47 

54 

59 

66 


(8.2832E - 09) 

(5.6681E - 09) 

(9.5469E - 09) 

(7.6937E - 09) 

4 

46 

53 

87 

65 


(9.8170E - 09) 

(8.4899E - 09) 

(7.9785E - 09) 

(1.6182E; -09) 

5 

43 

44 

56 

56 


(2.7604E - 09) 

(6.3987E - 09) 

(7.2821E - 09) 

(6.2222E; - 09) 

6 

36 

38 

105 

67 


(5.3518E - 09) 

(5.7851E - 09) 

(1.2954E - 10) 

(3.0156E - 10) 

7 

41 

43 

76 

77 


(3.7360E - 09) 

(3.1947E - 09) 

(1.1064E - 10) 

(4.9070E - 10) 

8 

38 

56 

86 

87 


(1.9864E - 09) 

(1.9460E -09) 

(4.8067E - 11) 

(4.6454E; - 10) 

9 

51 

51 

96 

59 


(4.6989E - 12) 

(2.2425E - 11) 

(2.8033E - 11) 

(6.7102E; - 09) 

10 

56 

56 

85 

65 


(5.5395E - 10) 

(1.0988E - 11) 

(2.7617E - 09) 

(1.3631E -09) 


Table 22 

Solitary wave generation of (I47p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of n to achieve a residual error below TOL = 10“®. The 
residual error (|2^ at the last computed iterate is in parenthesis; = 1, T = 0.999. 



U.C 70 I u.s u.a<i u.SH 

r r 


(a) (b) 

Eig. 23. Solitary wave generation of (1^71) . c* = 1. Number of iterations required to 
achieve a residual error (in Euclidean norm) below TOL = 10“® as function of T. 
(a) MPE, RRE, VEA and TEA. (b) AA-I and AA-IL 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

3 

53 

54 

57 

72 


(5.7252E - 09) 

(6.4717E - 09) 

(4.9757E - 09) 

(8.5993E - 09) 

4 

51 

47 

82 

65 


(6.5144E - 09) 

(6.5364E - 09) 

(1.2292E -09) 

(7.8352E - 10) 

5 

43 

49 

67 

68 


(1.4420E - 09) 

(5.7277E - 09) 

(9.3104E - 09) 

(3.3648E - 09) 

6 

50 

50 

105 

105 


(2.8215E - 12) 

(1.3107E - 10) 

(4.7779E - 12) 

(3.2453E - 09) 

7 

49 

50 

91 

122 


(8.1696E - 09) 

(5.7674E - 09) 

(1.4031E - 11) 

(6.5282E - 10) 

8 

55 

55 

103 

155 


(2.5994E - 10) 

(9.0192E - 09) 

(1.4896E - 10) 

(2.9075E - 10) 


Table 23 

Solitary wave generation of (|47p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error below TOL = 10“®. The 
residual error ([25|) at the last computed iterate is in parenthesis; c* = 1, T = 0.9999. 


nw 

AA-^nu)) 

AA-II(nri)) 

2 

23(2.5840E - 09) 

25(5.5612E - 09) 

3 


20(2.5548E - 09) 

4 

36(1.5755E - 09) 

18(8.3616E - 09) 

5 

18(6.8363E - 09) 

17(6.2358E - 09) 

6 

19(4.5057E - 09) 

17(6.4962E - 09) 

7 

19(3.3613E - 09) 

17(3.1330E - 09) 


Table 24 

Solitary wave generation of ()47p . Number of iterations required by AA-I and AA-II 
as function of nw to achieve a residual error below TOL = 10“®. The residual error 
(f25P at the last computed iterate is in parenthesis, c* = 1, T = 0.99. 


trated with the numerical generation of localized ground state solutions of the 
following generalized nonlinear Schrodinger equation 


iut + Uxx — V (x)u + \u\‘^u — 0.2|m|‘^m + u\ufu = 0, (53) 

with V{x) = —3.5sech^(a; + 1.5) — 3sech^(a: — 1.5) and u a real constant. Equa¬ 
tion fl53ll was studied in [81] (see also references therein) where a bifurcation 
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of solitary waves for z/ = z/c ~ 0.01247946 was analyzed. The bifurcation 
is of transcritical type with two tangentially connected branches of smooth 
solutions. This can be characterized by using the behaviour of the power 

/ CXD 

fi)dx, (54) 

-OO 

as function of p for any localized ground state solution u{x, t) = U (x, fi G 

M. The two branches are connected at some (po,-P(ho)) ~ (3.28,14.35). The 
numerical generation of localized ground state prohles of fl53|) with e-Petviashvili 
type methods was treated in [1] where the equation for the prohles U{x, /i) 


-fiU + u" - V{x)U + \U\^U - 0.2|17|^17 + v\U\^U = 0, 

was discretized by Fourier collocation techniques, leading to the system LUh = 
N{Uh) for the vector approximation Uh at the grid points Xj and where 


L = III - Dl + diag{V{xo),.. .,V{xm-i)), 

N{Uh) = N,{Uh) + miUh) + NsiUh) 

= {\Uh\.^) .Un - 0.2 (\Uh\.*) .Uh + K {\Uh\.^) .Uh. (55) 

The nonlinearity in fl5^ contains three homogeneous terms with degrees pi = 
3,P2 = 5,P3 = 7 and the e-Petviashvili method 






i=i 

_ ( {Lu,u) pj 

\{N{u),u)) ’ Pi-1’ 


J = 1,2,3, 


(56) 

(57) 


is applied. The iteration (l56|) . flFT)) will be considered as the method to be 
complemented with acceleration techniques. Finally, the quantity fl5T)) has 
been approximated by 


Pkiu,.) = hY.Ky 

j 


(58) 


The numerical illustration of this case takes p = 3.281, TOL = 10“^° and a 
superposition of squared hyperbolic secant functions as initial iteration. The 
numerical prohle Uh generated by fl56D . fl57)) is shown in Figure YIM a.). The 
corresponding value for (158|) is Ph{Uh) = 14.446162 and the poor performance 
of the method is made clear in Figure [24lf bi which displays the behaviour 
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Classical fixed point 

e-Petviashvili method ([56]), ([^ 

1.687048E-F00 

9.829607E-01 

9.834930E-01 

4.740069E-01 

4.793766E-01 

3.616157E-01 

3.747266E-01 

2.606251E-01-hl.734293E-01i 

1.979766E-01 

2.606251E-01-1.734293E-01i 

1.426054E-01 

1.764488E-01 


Table 25 

Six largest magnitude eigenvalues of the iteration matrices of classical fixed point 
algorithm and of e-Petviashvili method (I56p . (1571) for /r = 3.281 at the last com¬ 
puted iterate. The dominant egienvalue in the column on the right justifies the slow 
performance of the method. 



Fig. 24. Numerical ground state generation of (|53p with // = 3.281. (a) Approximate 
asymmetric profile; (b) Residual error as function of the number of iterations for 
the e-Petviashvili method (l56P . (f57P without acceleration. 

of the residual error fl2^ as function of the number of iterations and shows 
that the method requires n = 1023 iterations to achieve a residual error below 
TOL. (See Table 12^ first eigenvalue of the second column, to explain this slow 
behaviour.) 

The application of acceleration with VEM to this example is displayed in 
Table 1261 and Figure |25l The following points are emphasized: 

• The acceleration leads to a great improvement with respect to the e-Petviashvili 
method fISBp . flFTP . In order to have a residual error below TOL, the reduc¬ 
tion in the number of iterations is above 90%. 

• As in the previous examples, polynomial methods work better than e- 
algorithms. By comparing the two polynomial methods, MPE is more effi¬ 
cient. Its best performance requires a large number of k (which means a long 
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K 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

2 

185 

260 

1250 

187 


(5.3669E - 11) 

(9.9728E - 11) 

(9.5278E - 11) 

(5.5118E - 11) 

3 

135 

135 

209 

155 


(5.1486E - 11) 

(7.4460E - 11) 

(8.6758E - 11) 

(6.5522E - 11) 

4 

118 

64 

167 

109 


(6.3944E - 11) 

(7.6714E - 11) 

(8.7013E - 11) 

(5.6136E - 11) 

5 

64 

69 

75 

78 


(7.0759E - 11) 

(8.7100E - 11) 

(7.1288E - 11) 

(8.7351E - 11) 

6 

55 

65 

85 

795 


(8.4602E - 11) 

(4.5311E - 11) 

(7.6708E - 11) 

(6.9120E - 11) 

7 

53 

58 

80 

91 


(6.6155E - 11) 

(1.1617E - 11) 

(9.8698E - 11) 

(3.4911E - 11) 

8 

49 

58 

89 

70 


(5.6418E - 11) 

(9.0426E - 11) 

(9.0902E - 11) 

(6.9110E - 11) 

9 

52 

62 

82 

78 


(1.3763E - 11) 

(4.1423E - 11) 

(7.4598E - 11) 

(7.1601E - 11) 

10 

47 

67 

88 

106 


(6.1296E - 11) 

(4.6194E - 11) 

(6.9450E - 11) 

(4.9290E - 11) 


Table 26 

Ground state generation of (I53p . Number of iterations required by MPE, RRE, 
VEA and TEA as function of k to achieve a residual error below TOL = 

The residual error ([25|) at the last computed iterate is in parenthesis; fj. = 3.281. 
For the e-Petviashvili method (j56|] . (I57p without acceleration n = 1023 iterations 
are required for a residual error of 9.9939E — 11. 

cycle, above eight). On the other hand, the best results for the e-algorithms 
are obtained with moderate values of k, around hve. (This also happens in 
general in the previous examples.) 

• The value of /i considered for the experiments is close to the one corre¬ 
sponding to the bifurcation point, that is it is close to the tangential point 
of the two branches of solitary wave solutions. The computation of the quan¬ 
tity flSS]l for each acceleration, shown in Table |53 attempts to study the 
behaviour of the iterations close to the bifurcation. In most of the cases, 
the computed value coincides to that of the prohle generated by fl56H . (1571) 
without acceleration. In the case of MPE(8) and TEA(2)-TEA(6), the value 
changes to Ph{Uh) = 14.458882. This suggests that for these cases the accel- 
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K, 

MPE(k) 

RRE(k) 

VEA(k) 

TEA(k) 

2 

14.446162A -h 00 

14.446162E 00 

14.446162E -b 00 

14.458882E -b 00 

3 

14.446162A -h 00 

14.446162E -b 00 

14.446162E -b 00 

14.458882E -b 00 

4 

14.446162A 00 

14.446162E -b 00 

14.446162E -b 00 

14.458882E -b 00 

5 

14.446162A -h 00 

14.446162E -b 00 

14.446162E -b 00 

14.458882E -b 00 

6 

14.446162A -h 00 

14.446162E -b 00 

14.446162E -b 00 

14.458882E -b 00 

7 

14.446162A -h 00 

14.446162E -b 00 

14.446162E -b 00 

14.446162E -b 00 

8 

14.458882A -h 00 

14.446162E -b 00 

14.446162E -b 00 

14.446162E -b 00 

9 

14.446162A 00 

14.446162E -b 00 

14.446162E -b 00 

14.446162E -b 00 

10 

14.446162A -h 00 

14.446162E -b 00 

14.446162E -b 00 

14.446162E -b 00 


Table 27 


Ground state generation of (l53|) . Values of ([58l) for each iteration from Table [26l 
/i = 3.281. 

erated iteration converges to the profile of the upper branch while in most 
of the cases (including the one without acceleration) the limit prohle be¬ 
longs to the lower branch, |8T]. (Close to the bifurcation indeed, the form of 
the profiles is very similar, see Figures [261( a) and (b). Note however from 
Table 1301 that the dominant eigenvalue of the iteration matrix of (|56|1 . (157)1 
is above one. This and Table [25] may explain the convergence of this method 
to the profile with PhiJJh) = 14.4446162. 

• The comparison with the best choices of the VEM is illustrated in Figure 
l25l which compares the behaviour of the residual error (1251) as function of 
the number of iterations and of CPU time in seconds. The results reveal 
again the better performance of the polynomial methods when the residual 
error starts to be below 10“®. 

When the iteration fl56ll . fl57|) is accelerated with the AAM, we obtain the 
results displayed in Table [28j 

• The behaviour of the methods in this case looks similar to that of some 
previous examples as far as the general performance is concerned: they are 
competitive for moderate values of nw, with a better performance of AA-II, 
less affected by ill-conditioning. 

• In some cases the AMM approximate profiles (see Figures l26lf ci and (d)) 
which correspond to values of (1581) out of the branches. This uncertain be¬ 
haviour provides the main drawback of the methods. The spectral informa¬ 
tion for the two additional approximate profiles is given in Tables [29] and 
[371 The results suggest the lack of preservation of (1581) through the iterative 
process. 
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Fig. 25. Numerical ground state generation of (I53p with // = 3.281. Residual error 
(j25h as function of the number of iterations (a) and CPU time in seconds (b) for 
the e-Petviashvili method (f56]l . ([571) accelerated with MPE(IO), RRE(7), VEA(5) 
and TEA(5). 

5 Concluding remarks and future work 


In this paper we have studied numerically the use of acceleration techniques 
applied to hxed point algorithms of Petviashvili type to generate numerically 
traveling waves in nonlinear dispersive wave equations. The comparison has 
been established between vector extrapolation methods and Anderson accel¬ 
eration methods for different types of traveling waves. From the plethora of 
numerical experiments, our main conclusions are: 

• The use of acceleration techniques improves the performance of the Petvi- 
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nw 

AA-I(nt(;) 

P 

AA-II(nTi;) 

P 

1 

78(5.28027; - 11) 

14.4589 

59(5.4557.E - 11) 

14.4589 

2 

50(7.88877; - 11) 

14.4462 

51(5.5612.E - 09) 

14.4462 

3 

74(4.48147; - 11) 

14.4589 

28(1.7879.E - 11) 

3.9918 

4 

93(9.03117; - 11) 

14.4462 

32(2.6548.E - 11) 

3.9918 

5 

49(2.04217; - 11) 

3.9918 

64(1.2763.E - 11) 

3.9918 

6 

64(3.35407; - 11) 

3.9918 

80(3.0285.E - 11) 

3.9918 

7 

61(4.51367; - 11) 

14.4462 

29(9.1724.E - 11) 

3.9918 

8 

57(4.70257; - 11) 

14.4662 

54(3.7450.E - 12) 

9.7217 

9 



55(5.6993.E - 12) 

3.9918 

10 



108(3.4559.E - 11) 

3.9918 


Table 28 

Ground state generation of (l53]l with /r = 3.281. Number of iterations required by 
AA-I and AA-II as function of nw to achieve a residual error below TOL = 

The residual error (I25p at the last computed iterate is in parenthesis. 


Classical fixed point 

e-Petviashvili method (l56]l. (l57l) 

1.994420E-F00 

4.824994E-01 

4.828929E-01 

2.921950E-01-4.208258E-02i 

2.155227E-01 

2.921950E-01-h4.208258E-02i 

1.266822E-01 

1.270267E-01 

8.385958E-02 

8.457380E-02 

6.009441E-02 

6.014597E-02 


Table 29 

Six largest magnitude eigenvalues of the iteration matrices of classical fixed point 
algorithm and of e-Petviashvili method (f56l) . (f57|) for /r = 3.281 at the last computed 
iterate for P = 3.9918. 

ashvili type methods in all the cases. This improvement is observed in two 
main points: hrst, when the Petviashvili type method is convergent, the ac¬ 
celeration reduces the number of iterations in a relevant way. (In some cases, 
this is really important: in some one-dimensional problems, the reduction 
is at least of 50% and attains up to 75%.) On the other hand, the mech¬ 
anism of acceleration, especially in the case of VEM, allows to transform 
initially divergent sequences into convergent processes. This is particularly 
relevant in traveling waves with high oscillations. Furthermore, acceleration 
has been shown to be more efficient than other alternatives for some cases, 
like numerical continuation. 
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Classical fixed point 

e-Petviashvili method (f56]l. (fFr|) 

1.643665E-F00 

1.016836E-F00 

1.015912E-F00 

4.764502E-01 

4.862159E-01 

3.707178E-01 

3.756950E-01 

2.883210E-01-1.334802E-01i 

2.022740E-01 

2.883210E-01-hl.334802E-01i 

1.417784E-01 

1.748132E-01 


Table 30 

Six largest magnitude eigenvalues of the iteration matrices of classical fixed point 
algorithm and of e-Petviashvili method (|56|) . (|57|) for // = 3.281 at the last computed 
iterate for P = 14.4559. 


Classical fixed point 

e-Petviashvili method (f56]l. (IFr|) 

1.919934E-F00 

1.418879E-F00 

1.201956E-F00 

6.091959E-01-4.750746E-02i 

6.134143E-01 

6.091959E-01-h4.750746E-02i 

3.716346E-01 

4.517846E-01 

2.394216E-01 

2.884719E-01 

1.733253E-01 

1.799915E-01 


Table 31 

Six largest magnitude eigenvalues of the iteration matrices of classical fixed point 
algorithm and of e-Petviashvili method (f56l) . (f571) for /r = 3.281 at the last computed 
iterate for P = 9.7217. 

• In general, VEM provide better results and among them, polynomial meth¬ 
ods such as MPE and RRE are more efficient than e-algorithms like VEA 
and TEA, although in some convergent cases the acceleration in terms of 
the number of iterations is very similar among all the methods while in 
computational time the e-algorithms work better. 

• The AAM are competitive in some cases but they are mostly affected by 
ill-conditioning and a more computational effort due to their longer imple¬ 
mentation. The best results of these methods are obtained when generating 
numerically periodic traveling waves in some nonlinear dispersive systems 
and ground state prohles in NLS type equations, while their performance is 
poor when computing highly oscillatory traveling waves. Their application 
to the e-Petviashvili type methods suggests an uncertain behaviour with 
respect to relevant quantities of the problem through the iteration. 

The main question to that this comparative study has not been able to answer 
is in the authors’ opinion Ending a deeper understanding of the way how ac¬ 
celeration techniques (especially VEM) work on these problems. In particular. 
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(a) (b) 




(c) 


(d) 


Fig. 26. Numerical ground state generation of (I53p with ^ = 3.281. (a) Approxi¬ 
mate asymmetric profile; (a) P = 14.4462; (b) P = 3.9918; (c) P = 14.4559; (d) 
P = 9.7217. 

we miss some conclusions about the width of extrapolation (which is related 
to the extrapolation step for convergence) to be used a priori (if it is possible 
to do that). We have observed that this looks to be strongly dependent on the 
problem under study. This might be though a good starting point for a future 
research. 
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